Code to load tools and data
 # packages
    require(arm)
    require(data.table)
    require(effects)
    require(foreach)
    require(ggimage)
    require(ggplot2)
    require(ggpubr)
    require(ggsci)
    require(ggtext)
    require(grid)
    require(gtable)
    require(here)
    require(kableExtra)
    require(MASS)
    require(multcomp)
    require(optimx)
    require(performance)  
    require(PerformanceAnalytics)
    require(png)
    require(RColorBrewer)
    require(rmeta)
    require(rphylopic)
    require(scales)
    require(viridis)
 # constants
    save_plot = TRUE
    round_ = 3 # number of decimal places to round model coefficients
    nsim = 5000 # number of simulations to extract estimates and 95%CrI
    ax_lines = "grey60" # defines color of the axis lines
    #colors <- c("#999999", "#E69F00", "#56B4E9") #viridis(3)
    set.seed(42)
    #width_ = .7 # spacing between error bars
    #col_ = c(brewer.pal(n =12, name = "Paired"), 'grey30','grey80')
 # functions
    # to add images to panels
      annotation_custom2 <- function (grob, xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, data) {
        layer(data = data, stat = StatIdentity, position = PositionIdentity, 
            geom = ggplot2:::GeomCustomAnn,
            inherit.aes = TRUE, params = list(grob = grob, 
                                              xmin = xmin, xmax = xmax, 
                                              ymin = ymin, ymax = ymax))
      }
    # to remove ggplot components
        gtable_filter_remove <- function (x, name, trim = TRUE){
          matches <- !(x$layout$name %in% name)
          x$layout <- x$layout[matches, , drop = FALSE]
          x$grobs <- x$grobs[matches]
          if (trim) 
            x <- gtable_trim(x)
          x
        }
    # customized ggplot theme
        theme_MB = theme(  
                  title = element_text(size=8, colour="grey30"),
                  axis.line = element_blank(),
                  #axis.line = element_line(colour="grey70", size=0.25),
                  axis.title = element_text(size=7, colour="grey30"),
                  axis.title.y = element_text(vjust=3.5),
                  axis.title.x = element_text(vjust=1),
                  axis.text = element_text(size=6),#, vjust = 0.5, hjust=1),# margin=units(0.5,"mm")),
                  axis.ticks.length=unit(0.5,"mm"),
                  axis.ticks = element_line(colour = "grey70", size = 0.1),
                  #axis.ticks.margin,
                  
                  strip.text.x = element_text(size = 6, color="grey30",  margin=margin(1,1,1,1,"mm")), #grey50
                  strip.text.y = element_text(size = 6, color="grey30",  margin=margin(1,1,1,1,"mm")), #grey50
                  strip.background = element_rect(fill="grey99",colour="grey70", size=0.25),
                    #strip.background = element_blank(), 
                    #strip.text = element_blank(),
                  panel.spacing = unit(0, "mm"),
                  panel.background=element_blank(),
                  panel.border = element_rect(colour="grey70", size=0.1, fill = NA), #panel.border=element_blank(),
                  panel.grid = element_blank(),

                  legend.text=element_text(size=6),
                  legend.title=element_text(size=6),
                  legend.key = element_rect(colour = NA, fill = NA),
                  legend.key.height= unit(0.5,"line"),
                  legend.key.width = unit(0.25, "cm"),
                  legend.margin = margin(0,0,0,0, unit="cm"),
                  legend.box.margin = margin(l = -6), #legend.justification = c(-1,0),
                  legend.background = element_blank()
                  )  
    # for estimates
        est_out =function(model = m, label = "", nsim = 5000){
            bsim = sim(model, n.sim=5000) 
            v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
            ci = apply(bsim@fixef, 2, quantile, prob=c(0.025,0.975)) 
            sd = apply(bsim@fixef, 2, sd)
            data.table(predictor=rownames(coef(summary(model))),estimate=v, lwr=ci[1,], upr=ci[2,], sd = sd, model = paste(label, "N =", nobs(model)))
          }
    # change color
      change_col = function(replace_black, theimg) {
        r_b = col2rgb(replace_black) / 255
        #theimg[theimg == 1] <- 2
        for (i in 1:3) {
            theimg[,,i][theimg[,,i] == 0] <- r_b[i]
        }
        return(theimg)
        }
    # for Supplementary Table output based on sim
      m_out = function(model = m, type = "mixed", 
        name = "define", dep = "define", fam = 'Gaussian',
        round_ = 3, nsim = 5000, aic = FALSE, save_sim = here::here('Data/model_sim/'), back_tran = FALSE, perc_ = 1){
          # perc_ 1 = proportion or 100%
        bsim = sim(model, n.sim=nsim)  
        
        if(save_sim!=FALSE){save(bsim, file = paste0(save_sim, name,'.RData'))}
       
        if(type != "mixed"){
          v = apply(bsim@coef, 2, quantile, prob=c(0.5))
          ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 

          if(back_tran == TRUE & fam == "binomial"){
           v = perc_*plogis(v)
           ci = perc_*plogis(ci)
           }
          if(back_tran == TRUE & fam == "binomial_logExp"){
                v = perc_*(1-plogis(v))
                ci = perc_*(1-plogis(ci))
                ci = rbind(ci[2,],ci[1,])
               }

          if(back_tran == TRUE & fam == "Poisson"){
           v = exp(v)
           ci = exp(ci)
          }

         oi=data.frame(type='fixed',effect=rownames(coef(summary(model))),estimate=v, lwr=ci[1,], upr=ci[2,])
          rownames(oi) = NULL
          oi$estimate_r=round(oi$estimate,round_)
          oi$lwr_r=round(oi$lwr,round_)
          oi$upr_r=round(oi$upr,round_)
          if(perc_ == 100){
           oi$estimate_r = paste0(oi$estimate_r,"%")
           oi$lwr_r = paste0(oi$lwr_r,"%")
           oi$upr_r = paste0(oi$upr_r,"%")
          }
         x=data.table(oi[c('type',"effect", "estimate_r","lwr_r",'upr_r')]) 
       
        }else{
         v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
         ci = apply(bsim@fixef, 2, quantile, prob=c(0.025,0.975)) 

         if(back_tran == TRUE & fam == "binomial"){
          v = perc_*plogis(v)
          ci = perc_*plogis(ci)
         }
          if(back_tran == TRUE & fam == "binomial_logExp"){
                v = perc_*(1-plogis(v))
                ci = perc_*(1-plogis(ci))
                ci = rbind(ci[2,],ci[1,])
               }

          if(back_tran == TRUE & fam == "Poisson"){
            v = exp(v)
            ci = exp(ci)
         }

        oi=data.table(type='fixed',effect=rownames(coef(summary(model))),estimate=v, lwr=ci[1,], upr=ci[2,])
            rownames(oi) = NULL
            oi[,estimate_r := round(estimate,round_)]
            oi[,lwr_r := round(lwr,round_)]
            oi[,upr_r :=round(upr,round_)]
            if(perc_ == 100){
             oi[,estimate_r := paste0(estimate_r,"%")]
             oi[,lwr_r := paste0(lwr_r,"%")]
             oi[,upr_r := paste0(upr_r,"%")]
            }
         oii=oi[,c('type',"effect", "estimate_r","lwr_r",'upr_r')] 
        
         l=data.frame(summary(model)$varcor)
         l = l[is.na(l$var2),]
         l$var1 = ifelse(is.na(l$var1),"",l$var1)
         l$pred = paste(l$grp,l$var1)

         q050={}
         q025={}
         q975={}
         pred={}
         
         # variance of random effects
         for (ran in names(bsim@ranef)) {
           ran_type = l$var1[l$grp == ran]
           for(i in ran_type){
            q050=c(q050,quantile(apply(bsim@ranef[[ran]][,,ran_type], 1, var), prob=c(0.5)))
            q025=c(q025,quantile(apply(bsim@ranef[[ran]][,,ran_type], 1, var), prob=c(0.025)))
            q975=c(q975,quantile(apply(bsim@ranef[[ran]][,,ran_type], 1, var), prob=c(0.975)))
            pred= c(pred,paste(ran, i))
            }
           }
         # residual variance
         q050=c(q050,quantile(bsim@sigma^2, prob=c(0.5)))
         q025=c(q025,quantile(bsim@sigma^2, prob=c(0.025)))
         q975=c(q975,quantile(bsim@sigma^2, prob=c(0.975)))
         pred= c(pred,'Residual')

         ri=data.table(type='random %',effect=pred, estimate_r=round(100*q050/sum(q050)), lwr_r=round(100*q025/sum(q025)), upr_r=round(100*q975/sum(q975)))
           
         ri[lwr_r>upr_r, lwr_rt := upr_r]
         ri[lwr_r>upr_r, upr_rt := lwr_r]
         ri[!is.na(lwr_rt), lwr_r := lwr_rt]
         ri[!is.na(upr_rt), upr_r := upr_rt]
         ri$lwr_rt = ri$upr_rt = NULL

         ri[,estimate_r := paste0(estimate_r,'%')]
         ri[,lwr_r := paste0(lwr_r,'%')]
         ri[,upr_r := paste0(upr_r,'%')]
        
        x = data.table(rbind(oii,ri))
        }
        
        x[1, model := name]                                                                
        x[1, response := dep]                                                                
        x[1, error_structure := fam]      
        N = length(resid(model))                                                          
        x[1, N := N ]                                                                

        x=x[ , c('model', 'response', 'error_structure', 'N', 'type',"effect", "estimate_r","lwr_r",'upr_r')] 

        if (aic == TRUE){   
            x[1, AIC := AIC(update(model,REML = FALSE))] 
            }
        if (aic == "AICc"){
            aicc = AICc(model)
            x[1, AICc := aicc] 
        }
        if(type == "mixed" & nrow(x[type=='random %' & estimate_r =='0%'])==0){
          R2_mar = as.numeric(r2_nakagawa(model)$R2_marginal)
          R2_con = as.numeric(r2_nakagawa(model)$R2_conditional)
          x[1, R2_mar := R2_mar]
          x[1, R2_con := R2_con]
         }
        x[is.na(x)] = ""
        return(x)
      } 
    # model assumption function
      m_ass = function(name = 'define', mo = m0, dat = d, fixed = NULL, categ = NULL, trans = "none", spatial = TRUE, temporal = TRUE, PNG = TRUE, outdir = 'outdir', n_col=8, width_ = 10, height_ = 5.5){
       l=data.frame(summary(mo)$varcor)
       l = l[is.na(l$var2),]
       nt = if(temporal==TRUE){1}else{0}
       ns = if(spatial==TRUE){7}else{0}
       n = 3+nrow(l)-1+length(fixed)+length(categ) +  nt +  ns
     
       if(PNG == TRUE){
        png(paste(outdir,name, ".png", sep=""), width=width_,height=height_,units="in",res=600) # width = 6
        par(mfrow=c(4, n_col),tcl = -0.08, cex = 0.5, cex.main = 0.9,#ceiling(n/n_col),n_col)
            oma = c(1,1,2,1),
            mar = c(2, 2, 2, 1), mgp=c(1,0,0)
            )
         }else{
          dev.new(width=width_,height=height_)
          par(mfrow=c(4,n_col), tcl = -0.08, cex = 0.5, cex.main = 0.9,#ceiling(n/n_col),n_col)
            oma = c(1,1,2,1),
            mar = c(2, 2, 2, 1), mgp=c(1,0,0)
            )
        }
       
       scatter.smooth(fitted(mo),resid(mo),col='grey');abline(h=0, lty=2, col ='red')
       scatter.smooth(fitted(mo),sqrt(abs(resid(mo))), col='grey')
       qqnorm(resid(mo), main=list("Normal Q-Q Plot: residuals"),col='grey');qqline(resid(mo), col = 'red')
       #unique(l$grp[l$grp!="Residual"])
       for(i in unique(l$grp[l$grp!="Residual"])){
        #i = "mean_year"
        ll=ranef(mo)[names(ranef(mo))==i][[1]]
        if(ncol(ll)==1){
         qqnorm(ll[,1], main = paste(i,names(ll)[1]),col='grey',);qqline(ll[,1], col ='red')
         }else{
          qqnorm(ll[,1], main = paste(i,names(ll)[1]),col='grey');qqline(ll[,1], col ='red')
          qqnorm(ll[,2], main = paste(i,names(ll)[2]),col='grey');qqline(ll[,2], col ='red')
         }
        }
        
       # variables
         scatter={} 
         for (i in rownames(summary(mo)$coef)) {
              #i = "lat_abs"
            j=sub("\\).*", "", sub(".*\\(", "",i)) 
            scatter[length(scatter)+1]=j
          }
          x = data.frame(scatter=unique(scatter)[2:length(unique(scatter))],
                          log_ = grepl("log",rownames(summary(mo)$coef)[2:length(unique(scatter))]), stringsAsFactors = FALSE)
          for (i in 1:length(fixed)){
              jj =fixed[i]
              variable=dat[, ..jj][[1]]
              if(trans[i]=='log'){
              scatter.smooth(resid(mo)~log(variable),xlab=paste('log(',jj,')',sep=''), col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
              }else if(trans[i]=='abs'){
              scatter.smooth(resid(mo)~abs(variable),xlab=paste('abs(',jj,')',sep=''), col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
              }else if(trans[i]=='sin'){scatter.smooth(resid(mo)~sin(variable),xlab=paste('sin(',jj,')',sep=''), col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
              }else if(trans[i]=='cos'){scatter.smooth(resid(mo)~cos(variable),xlab=paste('cos(',jj,')',sep=''), col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
              }else{
              scatter.smooth(resid(mo)~variable,xlab=jj,col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
            }
           }
          
          if(length(categ)>0){
            for(i in categ){
               variable=dat[, ..i][[1]]
                boxplot(resid(mo)~variable, medcol='grey', whiskcol='grey', staplecol='grey', boxcol='grey', outcol='grey', xlab = i);abline(h=0, lty=3, lwd=1, col = 'red')
               }
          }     
              
        if(temporal == TRUE){
            acf(resid(mo), type="p", main=list("Temporal autocorrelation:\npartial series residual"))
            }
        if(spatial == TRUE){    
          spdata=data.frame(resid=resid(mo), x=dat$Lon, y=dat$Lat)
            spdata$col=ifelse(spdata$resid<0,rgb(83,95,124,100, maxColorValue = 255),ifelse(spdata$resid>0,rgb(253,184,19,100, maxColorValue = 255), 'red'))
            #cex_=c(1,2,3,3.5,4)
            cex_=c(1,1.5,2,2.5,3)
            spdata$cex=as.character(cut(abs(spdata$resid), 5, labels=cex_))
          plot(spdata$x, spdata$y,col=spdata$col, cex=as.numeric(spdata$cex), pch= 16, main=list('Spatial distribution of residuals', cex=0.8), xlab = 'Longitude', ylab = 'Latitude')
            legend("topright", pch=16, legend=c('>0','<0'), ,col=c(rgb(83,95,124,100, maxColorValue = 255),rgb(253,184,19,100, maxColorValue = 255)))
          plot(spdata$x[spdata$resid<0], spdata$y[spdata$resid<0],col=spdata$col[spdata$resid<0], cex=as.numeric(spdata$cex[spdata$resid<0]), pch= 16, main=list('residuals <0'), xlab = 'Longitude', ylab = 'Latitude')
          plot(spdata$x[spdata$resid>=0], spdata$y[spdata$resid>=0],col=spdata$col[spdata$resid>=0], cex=as.numeric(spdata$cex[spdata$resid>=0]), pch= 16, main=list('residual >=0'), xlab = 'Longitude', ylab = 'Latitude')

          # EU
          dat$res = resid(mo)
          spdata=data.frame(resid = dat$res[dat$Country!='Australia'], x=dat$Lon[dat$Country!='Australia'], y=dat$Lat[dat$Country!='Australia'])
          spdata$col=ifelse(spdata$resid<0,rgb(83,95,124,100, maxColorValue = 255),ifelse(spdata$resid>0,rgb(253,184,19,100, maxColorValue = 255), 'red'))
            #cex_=c(1,2,3,3.5,4)
            cex_=c(1,1.5,2,2.5,3)
            spdata$cex=as.character(cut(abs(spdata$resid), 5, labels=cex_))
          plot(spdata$x[spdata$resid<0], spdata$y[spdata$resid<0],col=spdata$col[spdata$resid<0], cex=as.numeric(spdata$cex[spdata$resid<0]), pch= 16, main=list('EU -  residuals <0'), xlab = 'Longitude', ylab = 'Latitude')
          plot(spdata$x[spdata$resid>=0], spdata$y[spdata$resid>=0],col=spdata$col[spdata$resid>=0], cex=as.numeric(spdata$cex[spdata$resid>=0]), pch= 16, main=list('EU residuals >=0)'), xlab = 'Longitude', ylab = 'Latitude')

          # Australia
          spdata=data.frame(resid = dat$res[dat$Country=='Australia'], x=dat$Lon[dat$Country=='Australia'], y=dat$Lat[dat$Country=='Australia'])
          spdata$col=ifelse(spdata$resid<0,rgb(83,95,124,100, maxColorValue = 255),ifelse(spdata$resid>0,rgb(253,184,19,100, maxColorValue = 255), 'red'))
            #cex_=c(1,2,3,3.5,4)
            cex_=c(1,1.5,2,2.5,3)
            spdata$cex=as.character(cut(abs(spdata$resid), 5, labels=cex_))
          plot(spdata$x[spdata$resid<0], spdata$y[spdata$resid<0],col=spdata$col[spdata$resid<0], cex=as.numeric(spdata$cex[spdata$resid<0]), pch= 16, main=list('Australia residuals <0'), xlab = 'Longitude', ylab = 'Latitude')
          plot(spdata$x[spdata$resid>=0], spdata$y[spdata$resid>=0],col=spdata$col[spdata$resid>=0], cex=as.numeric(spdata$cex[spdata$resid>=0]), pch= 16, main=list('Australia residuals >=0'), xlab = 'Longitude', ylab = 'Latitude')
          }
       
       mtext(stringr::str_wrap(paste(paste0(name," model: "), slot(mo,"call")[1],'(',slot(mo,"call")[2],sep=''), width = ceiling(nchar(paste(slot(mo,"call")[1],'(',slot(mo,"call")[2],sep=''))/2)+10), side = 3, line = 0, cex=0.5,outer = TRUE, col = 'darkblue') #ceiling(nchar(paste(slot(mo,"call")[1],'(',slot(mo,"call")[2],sep=''))/2)
       if(PNG==TRUE){dev.off()}
      }
    
 # data
    ph  =  fread(here::here('Data/phylopic.txt'))
    setnames(ph, old = c('Name', 'Code'), new = c('genus2', 'uid'))

    t = fread(here::here('Data/taxonomy.txt'))
    
    g = fread(here::here('Data/google_mobility.txt')) #fwrite(d, here::here('Data/data.txt'), sep ='\t')
    g[, Year := as.integer(substring(date, nchar(date)-3, nchar(date)))]
    g[nchar(date)==9, date:=paste0('0',date)]
    g[, date_ :=as.Date(date, format = '%d.%m.%Y')]
    g[, Day :=yday(date_)]
    g[country_region!='Australia', Day := Day-92 +1] # 1 April = start of breeding season (1st day) = 92 day of the year 
    g[country_region=='Australia', Day := Day-228 +1] # 15 Augusst = start of breeding season (1st day) = 228 day of the year 
    setnames(g, old = 'country_region', new ='Country')

    d = fread(here::here('Data/data_corrected.txt')) #fwrite(d, here::here('Data/data.txt'), sep ='\t')
    # add data and weekdays
    x = fread(here::here("Data/date.txt"))[, .(IDObs, Date_corr)]
    x[, date_:=as.Date(Date_corr, "%d.%m.%Y" )]
    x[, weekday := weekdays(date_)]
    d =  merge(d, x[, .(IDObs, date_, weekday)], by = "IDObs")
    #d[, date_ := as.Date(Day, origin = paste0(Year, '-01-01'))]

    # adjust correct assignment of season (Year) for Australia
    d[Country == 'Australia' & Year == 2020 & Covid == 0, Year:=2019]
    d[Country == 'Australia' & Year == 2021 & Day>139, Year:=2020]
    d[Country == 'Australia' & Year == 2022 & Day>139, Year:=2021]

    d = d[order(Year, IDLocality, Day, Hour)]
    d[Country %in% c("Czech_Republic", "Czech Republic"), Country := "Czechia"]
    d[, genus := sub("_.*", "", Species)]
    d[, sp_day_year := paste(Year, Species, Day, sep="_")]
    d[, sp_loc := paste(Species, IDLocality, sep="_")]
    d[, sp_country := paste(Species, Country, sep="\n")]
    d[, rad:=(2*pi*Hour) / 24]
    d[, Day_:= Day]
    #d[, FID_z := scale(FID), by = Species]
    #d[, SD_z := scale(SD), by = Species]
    d[, FID_ln := log(FID)]
    d[, SD_ln := log(SD)]
    d[, body_ln := log(BodyMass)]
    d[, flock_ln := log(FlockSize)]
    d[, weekday := weekdays(date_)]
    #d[Country == 'Australia', Day_:= abs(Day - 189)]
    
    d1 = d[Covid == 1, .N, by = Species]
    d2 = d[Covid == 0, .N, by = Species]
    setnames(d1, old = 'N', new ='N_during')
    setnames(d2, old = 'N', new ='N_before')
    dd = merge(d1,d2) # species with data before and during
    da = merge(d1,d2, all = TRUE)

    d1p = d[Country!='Poland' & Covid == 1, .N, by = Species]
    d2p = d[Country!='Poland' & Covid == 0, .N, by = Species]
    setnames(d1p, old = 'N', new ='N_during')
    setnames(d2p, old = 'N', new ='N_before')
    ddp = merge(d1p,d2p) # species with data before and during, but without Poland data

    p1 = d[Covid == 1, .N, by = .(IDLocality, Species)]
    p2 = d[Covid == 0, .N, by = .(IDLocality, Species)]
    setnames(p1, old = 'N', new ='N_during')
    setnames(p2, old = 'N', new ='N_before')
    pp = merge(p1,p2)  # species-localities with data before and during
    pa = merge(p1,p2, all = TRUE)

    p1p = d[Country!='Poland' & Covid == 1, .N, by = .(IDLocality, Species)]
    p2p = d[Country!='Poland' & Covid == 0, .N, by = .(IDLocality, Species)]
    setnames(p1p, old = 'N', new ='N_during')
    setnames(p2p, old = 'N', new ='N_before')
    ppp = merge(p1p,p2p)  # species-localities with data before and during,but without Poland data 
   
    p1p4 = d[Year!=2014 & Country!='Poland' & Covid == 1, .N, by = .(IDLocality, Species)]
    p2p4 = d[Year!=2014 &Country!='Poland' & Covid == 0, .N, by = .(IDLocality, Species)]
    setnames(p1p4, old = 'N', new ='N_during')
    setnames(p2p4, old = 'N', new ='N_before')
    ppp4 = merge(p1p4,p2p4) # species-localities with data before and during, but without Poland & 2014 data

    s = d[Covid == 1]
    s[, Nsp := .N, by ='Species']
    s[, sp := gsub('[_]', ' ', Species)]
    # add google mobility
    s = merge(s, g[,.(Country,  date_, parks_percent_change_from_baseline)], all.x = TRUE)
    s[, year_ := as.character(Year)]  

    ss = s[!is.na(parks_percent_change_from_baseline)]
    ss[, country_year := paste(Country, Year)] #table(paste(s$Country, s$Year))   
    ss[parks_percent_change_from_baseline<0, google := 'before_zero']
    ss[parks_percent_change_from_baseline>0, google := 'after_zero']
    ss[, sp_country_google:= paste(sp_country, google)]

    g[, weekday := weekdays(date_)]

correlations among predictors

  d[, sin_rad:=sin(rad)]
  d[, cos_rad := cos(rad)]

  dp = d[, c("SD_ln", "flock_ln", "body_ln", "sin_rad", "cos_rad", "Temp", "Day")]
  setnames(dp, old = c("SD_ln", "flock_ln", "body_ln", "sin_rad", "cos_rad", "Temp", "Day"), new = c("Starting distance\nln(m)", "Flock size\nln(m)", "Body mass\nln(m)", "Sine\nof radians", "Cosine\nof radians", "Temperature\n°C", "Day"))

  chart.Correlation(dp, histogram = TRUE, pch = 19, alpha = 0.5)
  mtext("Single observations", side = 3, line = 3)

  if(save_plot==TRUE){
  png(here::here("Outputs/Fig_S1_rev.png"), width = 19, height = 19, units = "cm", bg = "transparent", res = 600)
  chart.Correlation(dp, histogram = TRUE, pch = 19, alpha = 0.5)
  mtext("Single observations", side = 3, line = 3)
  dev.off()
  }
## quartz_off_screen 
##                 2

prediction for FID ~ Period

# predictions
# full
ms <- lmer(scale(log(FID)) ~
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1| genus) + (1| Species) + (1 | sp_day_year) + (scale(Covid) | Country) + (scale(Covid) | IDLocality) + (1 | sp_loc),
data = d, REML = FALSE,
control <- lmerControl(
    optimizer = "optimx", optCtrl = list(method = "nlminb")
)
)
est_ms <- est_out(ms, "All: (1|Year) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (1|sp_loc)")
est_ms[, control_for_starting_distance := "yes"]

mx <- lmer(scale(log(FID)) ~
    #scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | Country) + (scale(Covid) | IDLocality) + (1 | sp_loc),
    data = d, REML = FALSE,
        control = lmerControl( 
            optimizer ='optimx', optCtrl=list(method='nlminb')) 
        )  
# (Covid|IDLocality) +
est_mx <- est_out(mx, "All: (1|Year) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (1|sp_loc)")
est_mx[, control_for_starting_distance := "no"]

# CZ - singular fits only due to genera estimated as zero (removing it changes no results)
  cs <- lmer(scale(log(FID)) ~
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (scale(Covid)| IDLocality) + (1|sp_loc),
    #(1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
    data = d[Country == "Czechia"], REML = FALSE
    )
  est_cs <- est_out(cs, "Czechia: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
  est_cs[, control_for_starting_distance := "yes"]
  
  cx <- lmer(scale(log(FID)) ~
    #scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (scale(Covid)| IDLocality) + (1|sp_loc),
    #(1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
    data = d[Country == "Czechia"], REML = FALSE
    )
  est_cx <- est_out(cx, "Czechia: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
  est_cx[, control_for_starting_distance := "no"]

# FI
fs <- lmer(scale(log(FID)) ~
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (scale(Covid)| IDLocality) + (1|sp_loc),
data = d[Country == "Finland"], REML = FALSE
)
est_fs <- est_out(fs, "Finland: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_fs[, control_for_starting_distance := "yes"]

fx <- lmer(scale(log(FID)) ~
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
     (1 | Year) + (1 | weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (scale(Covid)| IDLocality) + (1|sp_loc),
data = d[Country == "Finland"], REML = FALSE
)
est_fx <- est_out(fx, "Finland: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_fx[, control_for_starting_distance := "no"]

# HU - singular fits only due to sp_loc estimated as zero (removing it changes no results)
hs <- lmer(scale(log(FID)) ~
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Hungary"], REML = FALSE
)
est_hs <- est_out(hs, "Hungary: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_hs[, control_for_starting_distance := "yes"]

hx <- lmer(scale(log(FID)) ~
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Hungary"], REML = FALSE
)
est_hx <- est_out(hx, "Hungary: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_hx[, control_for_starting_distance := "no"]

# AU - singular fits only due to Year and random slope estimated as zero (removing those changes no results)
as <- lmer(scale(log(FID)) ~
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
     (1 | Year) +(1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Australia"], REML = FALSE
)
est_as <- est_out(as, "Australia: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_as[, control_for_starting_distance := "yes"]

ax <- lmer(scale(log(FID)) ~
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Australia"], , REML = FALSE,
        control = lmerControl( 
            optimizer ='optimx', optCtrl=list(method='nlminb')) 
)
est_ax <- est_out(ax, "Australia: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_ax[, control_for_starting_distance := "no"]

# PL 
ps <- lmer(scale(log(FID)) ~
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Poland"], REML = FALSE
)
est_ps <- est_out(ps, "Poland: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)")
est_ps[, control_for_starting_distance := "yes"]

px <- lmer(scale(log(FID)) ~
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Poland"], REML = FALSE
)
est_px <- est_out(px, "Poland: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)")
est_px[, control_for_starting_distance := "no"]

  # combine
    est_ms[, Country := 'All\n(mixed model)']
    est_mx[, Country := "All\n(mixed model)"]
    est_as[, Country := "Australia"]
    est_ax[, Country := "Australia"]
    est_cs[, Country := "Czechia"]
    est_cx[, Country := "Czechia"]
    est_hs[, Country := "Hungary"]
    est_hx[, Country := "Hungary"]
    est_ps[, Country := "Poland"]
    est_px[, Country := "Poland"]
    est_fs[, Country := "Finland"]
    est_fx[, Country := "Finland"]

    o = rbind(est_ms, est_mx, 
            est_as, est_ax, 
            est_cs, est_cx, 
            est_hs, est_hx,
            est_ps, est_px, 
            est_fs, est_fx)
    save(o, file = here::here('Data/dat_est_rev.Rdata'))
 
  #AIC(hs, hx)
  #AIC(cs, cx)
  #0AIC(ps, px)
  #AIC(as, ax)
  #AIC(fs, fx)

# supplementary Sable S_fid_period  
  ms_out <- m_out(name = "Table S_FID_country - full a", dep = "Escape distance", model = ms, nsim = 5000)
  mx_out <- m_out(name = "Table S_FID_country - full b", dep = "Escape distance", model = mx, nsim = 5000)
  cs_out <- m_out(name = "Table S_FID_country - CZ a", dep = "Escape distancey", model = cs, nsim = 5000)
  cx_out <- m_out(name = "Table S_FID_country - CZ b", dep = "Escape distance", model = cx, nsim = 5000)
  fs_out <- m_out(name = "Table S_FID_country - FI a", dep = "Escape distance", model = fs, nsim = 5000)
  fx_out <- m_out(name = "Table S_FID_country - FI b", dep = "Escape distance", model = fx, nsim = 5000)
  hs_out <- m_out(name = "Table S_FID_country - HU a", dep = "Escape distance", model = hs, nsim = 5000)
  hx_out <- m_out(name = "Table S_FID_country - HU b", dep = "Escape distance", model = hx, nsim = 5000)
  as_out <- m_out(name = "Table S_FID_country - AU a", dep = "Escape distance", model = as, nsim = 5000)
  ax_out <- m_out(name = "Table S_FID_country - AU b", dep = "Escape distance", model = ax, nsim = 5000)
  ps_out <- m_out(name = "Table S_FID_country - PL a", dep = "Escape distance", model = ps, nsim = 5000)
  px_out <- m_out(name = "Table S_FID_country - PL b", dep = "Escape distance", model = px, nsim = 5000)
  
  out_FID_c <- rbind(fs_out, fx_out, ps_out, px_out, cs_out, cx_out, hs_out, hx_out, as_out, ax_out, fill = TRUE)
  out_FID_c[is.na(out_FID_c)] <- ""
  out_FID_c[, effect := gsub("scale\\(Covid\\)", "Period", effect)]
  out_FID_c[, effect := gsub("Year", "year", effect)]
  out_FID_c[, effect := gsub("scale\\(SD\\)", "starting distance", effect)]
  out_FID_c[, effect := gsub("Temp", "temperaturre", effect)]
  out_FID_c[, effect := gsub("FlockSize", "flock size", effect)]
  out_FID_c[, effect := gsub("BodyMass", "body mass", effect)]
  out_FID_c[, effect := gsub("Species", "species", effect)]
  out_FID_c[, effect := gsub("sp_day_year", "species within day & year", effect)]
  out_FID_c[, effect := gsub("IDLocality", "site", effect)]
  out_FID_c[, effect := gsub("sp_loc", "species within site", effect)]
  fwrite(file = here::here("Outputs/Table_S_FID_c.csv"), out_FID_c)
load(here::here("Data/dat_est_rev.Rdata"))
o[predictor %in% c("scale(Covid)"), predictor := "Period"]
oo <- o[predictor %in% c("Period")]
oo[, N:=as.numeric(sub('.*N = ', '', model))]
# add meta-analytical mean
  oo_s = oo[control_for_starting_distance == 'yes']
  met = summary(meta.summaries(d = oo_s$estimate, se = oo_s$sd, method = "fixed", weights = oo_s$N))$summci
  oo_met = data.table(predictor = "Period", estimate = met[2], lwr = met[1], upr = met[3], sd = NA, model = NA, control_for_starting_distance = "yes", Country = "Combined\n(metanalytical)", N = NA)

  oo_sx = oo[control_for_starting_distance == "no"]
  metx = summary(meta.summaries(d = oo_sx$estimate, se = oo_sx$sd, method = "fixed", weights = oo_sx$N))$summci
  oo_metx = data.table(predictor = "Period", estimate = metx[2], lwr = metx[1], upr = metx[3], sd = NA, model = NA, control_for_starting_distance = "no", Country = "Combined\n(metanalytical)", N = NA)
  
  oo = rbind(oo, oo_met, oo_metx)
    
oo[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia", "Combined\n(metanalytical)", "All\n(mixed model)")))]

width_ <- .5 # spacing between error bars

#col_ <- c(brewer.pal(n = 12, name = "Paired"), "grey30", "grey80")
#Tol_bright <- c("#EE6677", "#228833", "#4477AA", "#CCBB44", "#66CCEE", "#AA3377", "#BBBBBB")
#Tol_muted <- c("#88CCEE", "#44AA99", "#117733", "#332288", "#DDCC77", "#999933", "#CC6677", "#882255", "#AA4499", "#DDDDDD")
#Tol_light <- c("#BBCC33", "#AAAA00", "#77AADD", "#EE8866", "#EEDD88", "#FFAABB", "#99DDFF", "#44BB99", "#DDDDDD")

# From Color Universal Design (CUD): https://jfly.uni-koeln.de/color/
#Okabe_Ito <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#000000")
#col_ = Okabe_Ito[7:1]
# JAMA and LocusZoom modified order
#col_ =  c("#374E55FF", "#374E55FF", "#DF8F44FF", "#79AF97FF", "#00A1D5FF", "#B24745FF",  "#80796BFF") #"#6A6599FF",
#col_ <- c("#357EBDFF", "#9632B8FF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#D43F3AFF", "#D43F3AFF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
col_ = c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
#show_col(col_)

#p = 
ggplot(oo, aes(x = estimate, y = Country, col = Country, shape = control_for_starting_distance)) +
    geom_vline(xintercept = 0, color = "grey", linetype = "dotted") +
    geom_errorbarh(aes(xmin = lwr, xmax = upr), height = 0, position = ggstance::position_dodgev(width_)) +
    # geom_point(position = ggstance::position_dodgev(.6)) +
    geom_point(position = position_dodge(width = width_), bg = "white", size = 1.1) +
    # scale_color_viridis(discrete=TRUE, begin=0, end = 0.5)  +
    # scale_fill_viridis(discrete=TRUE, begin=0, end = 0.5) +
    # geom_text( aes(x = n_pos,label = N), vjust = 0, size = 1.75, position = ggstance::position_dodgev(width_))+ # 3 positions for 3 bars
    # annotate("text", x=log10(3), y=85, label= "Used", col = "grey30", size = 2.5)+

    scale_shape_manual(name = "Controlled for\nstarting distance", guide = guide_legend(reverse = TRUE), values = c(21, 19)) +
    #scale_color_jama(guide = "none")+ #, palette = 'light'
    scale_color_manual(guide = "none", values = col_) + #guide_legend(reverse = TRUE)
    scale_x_continuous(breaks = round(seq(-0.6, 1.2, by = 0.3), 1)) +
    ylab("") +
    xlab("Standardized effect size of period\n[on flight initiation distance]") +
    # coord_cartesian(xlim = c(-.15, .15)) +
    # scale_x_continuous(breaks = round(seq(-.15, .15, by = 0.05),2)) +
    theme_bw() +
    theme(
        legend.position = "right",
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 6),
        # legend.spacing.y = unit(0.1, 'cm'),
        legend.key.height = unit(0.5, "line"),
        legend.margin = margin(0, 0, 0, 0),
        # legend.position=c(0.5,1.6),
        plot.title = element_text(color = "grey", size = 7),
        plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 0.5, unit = "pt"),
        panel.grid = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = ax_lines, size = 0.25),
        axis.line.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
        # axis.text.x = element_text()
        axis.ticks.length = unit(1, "pt"),
        axis.text.x = element_text(, size = 6),
        axis.text.y = element_text(colour = "black", size = 7),
        axis.title = element_text(size = 7)
    )

ggsave(here::here("Outputs/Fig_rev_width_CustomLocusZoom_v2.png"), width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_CustomJAMAv1.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_Okabe_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_UChicago_v3.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_d3_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_nejm_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_jama_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_jco_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_npg_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)

# show used colors
#gg <- ggplot_build(p)
#col_ = unique(gg$data[[3]]["colour"])$colour
#show_col(col_)

Fig. 1 | The effect size with 95%CIs of Period (0 – before, 1 – during shutdowns) on flight initiation distance (ln-transformed). The general model on all data was controlled for starting distance (ln-transformed; filled circles) or not (empty circles), flock size (ln-transformed), temperature (also a proxy for a day within the breeding season: r Pearson = 0.48; Fig. S1) and time of day. To account for circular properties of time, time was transformed into radians (2 × time × π/24) and fitted as sine and cosine of radians (Bulla et al. 2016). All continuous variables were standardised by subtracting the mean and dividing by the standard deviation. The multicollinearity was small as correlations between predictors were weak (r<XX, Fig. S1). To account for the non-independence of data points (Schielzeth & Forstmeier 2009; Barr et al. 2013), we attempted to fit random intercepts of year, weekday, genus, species, species at a given day and year, country, site, and species within a site, while fitting Period as random slope within site. Note that fitting Period as random slope at other random intercepts produces similar results (see Fig Sxx). We used this approach with a full dataset with all observations (n = 6369), as well as with conservative datasets, one with at least five observations per species and Period (i.e. at least five observations before and five during the COVID-19 shutdowns; n = 5260), the other with at least 10 observations per species and each Period (n = 5106). In other words, conservative datasets included species sampled in the both Periods. Albeit random structure of the full model accounts for the potential country specific effect (country accounted for xx% of variance, Table xx), depicted are also estimates from country specific models, with same predictors and random structure as the full model (excluding the random intercept of country and in case of Poland also site, as specific sites were not noted in Poland). Using meta.summaries function from rmeta R-package we used the country estimates, their standard deviation, and sample size per country to estimate the meta-analytical mean, which reflects the results based on the full mixed effect model. Also, albeit the response of humans to shutdowns were likely country specific and presence of humans in parks might have increased in some countries, decreased in others or did not change, the escape distances do not reflect such changes. In other words, the birds are either inflexible or the change in human behavior due to shutdowns was not strong enought, which might have been the case - see Fig. ZZ - decide what to show in this figure (I think 2020-2022 changes in Google Mobility for each country including the fits and refering to supplementary figure)

Table FID_c | Escape distance in relations to Period, given country

out_FID_c$R2_mar = out_FID_c$R2_con = NULL
out_FID_c %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
model response error_structure N type effect estimate_r lwr_r upr_r
Table S_FID_country - FI a Escape distance Gaussian 1019 fixed (Intercept) -0.035 -0.321 0.231
fixed scale(log(SD)) 0.282 0.228 0.336
fixed scale(log(flock size)) -0.049 -0.101 0.001
fixed scale(log(body mass)) 0.091 -0.066 0.24
fixed scale(sin(rad)) 0.085 0.022 0.147
fixed scale(cos(rad)) 0.031 -0.028 0.088
fixed scale(temperaturre) -0.079 -0.148 -0.009
fixed Period -0.15 -0.337 0.031
random % species within day & year (Intercept) 10% 9% 11%
random % species within site (Intercept) 6% 5% 7%
random % site (Intercept) 1% 0% 3%
random % site Period 1% 0% 3%
random % species (Intercept) 4% 3% 5%
random % genus (Intercept) 10% 7% 14%
random % weekday (Intercept) 0% 0% 1%
random % year (Intercept) 5% 1% 11%
random % Residual 61% 50% 72%
Table S_FID_country - FI b Escape distance Gaussian 1019 fixed (Intercept) -0.125 -0.468 0.221
fixed scale(log(flock size)) -0.018 -0.07 0.034
fixed scale(log(body mass)) 0.2 0.03 0.375
fixed scale(sin(rad)) 0.093 0.026 0.161
fixed scale(cos(rad)) 0.029 -0.032 0.089
fixed scale(temperaturre) -0.076 -0.148 -0.005
fixed Period -0.132 -0.366 0.101
random % species within day & year (Intercept) 9% 7% 10%
random % species within site (Intercept) 7% 6% 7%
random % site (Intercept) 1% 1% 2%
random % site Period 1% 1% 2%
random % species (Intercept) 7% 5% 8%
random % genus (Intercept) 10% 7% 13%
random % weekday (Intercept) 1% 0% 1%
random % year (Intercept) 9% 2% 16%
random % Residual 56% 44% 68%
Table S_FID_country - PL a Escape distance Gaussian 762 fixed (Intercept) 0.145 -0.115 0.404
fixed scale(log(SD)) 0.479 0.419 0.538
fixed scale(log(flock size)) 0.012 -0.042 0.066
fixed scale(log(body mass)) -0.07 -0.172 0.031
fixed scale(sin(rad)) -0.035 -0.117 0.05
fixed scale(cos(rad)) -0.012 -0.101 0.078
fixed scale(temperaturre) -0.111 -0.182 -0.04
fixed Period 0.119 -0.112 0.353
random % species within day & year (Intercept) 8% 8% 8%
random % species (Intercept) 15% 12% 17%
random % genus (Intercept) 0% 0% 0%
random % weekday (Intercept) 0% 0% 1%
random % year (Intercept) 11% 7% 17%
random % Residual 66% 57% 72%
Table S_FID_country - PL b Escape distance Gaussian 762 fixed (Intercept) 0.246 -0.075 0.558
fixed scale(log(flock size)) 0.054 -0.009 0.117
fixed scale(log(body mass)) 0.046 -0.097 0.19
fixed scale(sin(rad)) -0.014 -0.112 0.077
fixed scale(cos(rad)) -0.001 -0.104 0.099
fixed scale(temperaturre) -0.104 -0.187 -0.02
fixed Period 0.245 0.001 0.498
random % species within day & year (Intercept) 11% 10% 11%
random % species (Intercept) 26% 22% 29%
random % genus (Intercept) 1% 1% 2%
random % weekday (Intercept) 0% 0% 0%
random % year (Intercept) 8% 5% 14%
random % Residual 53% 46% 60%
Table S_FID_country - CZ a Escape distancey Gaussian 2013 fixed (Intercept) 0.111 -0.257 0.481
fixed scale(log(SD)) 0.392 0.351 0.437
fixed scale(log(flock size)) 0.005 -0.03 0.037
fixed scale(log(body mass)) 0.17 0.038 0.308
fixed scale(sin(rad)) 0.05 0.002 0.099
fixed scale(cos(rad)) 0.054 0.009 0.097
fixed scale(temperaturre) 0.043 -0.008 0.093
fixed Period 0.036 -0.294 0.361
random % species within day & year (Intercept) 3% 2% 3%
random % species within site (Intercept) 3% 2% 3%
random % species (Intercept) 13% 12% 13%
random % site (Intercept) 2% 0% 3%
random % site Period 2% 0% 3%
random % genus (Intercept) 12% 11% 11%
random % weekday (Intercept) 0% 0% 0%
random % year (Intercept) 15% 1% 33%
random % Residual 52% 34% 70%
Table S_FID_country - CZ b Escape distance Gaussian 2013 fixed (Intercept) -0.043 -0.969 0.867
fixed scale(log(flock size)) 0.005 -0.031 0.04
fixed scale(log(body mass)) 0.231 0.078 0.382
fixed scale(sin(rad)) 0.14 0.086 0.195
fixed scale(cos(rad)) 0.123 0.075 0.172
fixed scale(temperaturre) 0.035 -0.027 0.098
fixed Period 0.223 -0.642 1.1
random % species within day & year (Intercept) 2% 1% 5%
random % species within site (Intercept) 1% 1% 2%
random % species (Intercept) 6% 4% 10%
random % site (Intercept) 1% -1% 1%
random % site Period 1% -1% 1%
random % genus (Intercept) 8% 5% 11%
random % weekday (Intercept) 1% 0% 1%
random % year (Intercept) 54% 17% 73%
random % Residual 26% 13% 57%
Table S_FID_country - HU a Escape distance Gaussian 1456 fixed (Intercept) 0.097 -0.174 0.384
fixed scale(log(SD)) 0.485 0.438 0.531
fixed scale(log(flock size)) -0.012 -0.051 0.031
fixed scale(log(body mass)) -0.065 -0.184 0.059
fixed scale(sin(rad)) -0.077 -0.133 -0.019
fixed scale(cos(rad)) -0.011 -0.06 0.037
fixed scale(temperaturre) -0.023 -0.083 0.037
fixed Period -0.091 -0.271 0.097
random % species within day & year (Intercept) 5% 4% 5%
random % species within site (Intercept) 0% 0% 0%
random % species (Intercept) 3% 2% 3%
random % genus (Intercept) 7% 5% 8%
random % site (Intercept) 2% 0% 8%
random % site Period 2% 0% 8%
random % weekday (Intercept) 1% 0% 1%
random % year (Intercept) 5% 2% 8%
random % Residual 76% 59% 84%
Table S_FID_country - HU b Escape distance Gaussian 1456 fixed (Intercept) 0.078 -0.336 0.494
fixed scale(log(flock size)) -0.011 -0.057 0.035
fixed scale(log(body mass)) 0.187 0.006 0.376
fixed scale(sin(rad)) -0.111 -0.18 -0.044
fixed scale(cos(rad)) -0.007 -0.06 0.05
fixed scale(temperaturre) -0.025 -0.095 0.043
fixed Period -0.105 -0.307 0.109
random % species within day & year (Intercept) 7% 5% 8%
random % species within site (Intercept) 0% 0% 0%
random % species (Intercept) 2% 2% 2%
random % genus (Intercept) 18% 13% 18%
random % site (Intercept) 1% 0% 12%
random % site Period 1% 0% 12%
random % weekday (Intercept) 1% 0% 1%
random % year (Intercept) 5% 3% 7%
random % Residual 63% 42% 74%
Table S_FID_country - AU a Escape distance Gaussian 1119 fixed (Intercept) -0.064 -0.202 0.081
fixed scale(log(SD)) 0.502 0.452 0.553
fixed scale(log(flock size)) 0.02 -0.026 0.068
fixed scale(log(body mass)) -0.066 -0.163 0.029
fixed scale(sin(rad)) -0.032 -0.099 0.036
fixed scale(cos(rad)) -0.026 -0.093 0.041
fixed scale(temperaturre) 0.009 -0.045 0.063
fixed Period -0.009 -0.08 0.062
random % species within day & year (Intercept) 7% 6% 7%
random % species within site (Intercept) 13% 11% 13%
random % species (Intercept) 14% 11% 14%
random % genus (Intercept) 2% 1% 2%
random % site (Intercept) 0% -1% 8%
random % site Period 0% -1% 8%
random % weekday (Intercept) 0% 0% 0%
random % year (Intercept) 0% 0% 0%
random % Residual 64% 51% 68%
Table S_FID_country - AU b Escape distance Gaussian 1119 fixed (Intercept) -0.18 -0.432 0.058
fixed scale(log(flock size)) 0.051 -0.002 0.104
fixed scale(log(body mass)) 0.055 -0.073 0.184
fixed scale(sin(rad)) -0.017 -0.098 0.063
fixed scale(cos(rad)) -0.006 -0.08 0.069
fixed scale(temperaturre) 0.004 -0.058 0.067
fixed Period 0.156 -0.007 0.331
random % species within day & year (Intercept) 8% 5% 10%
random % species within site (Intercept) 12% 8% 15%
random % species (Intercept) 17% 14% 17%
random % genus (Intercept) 5% 4% 5%
random % site (Intercept) 1% -10% 15%
random % site Period 1% -10% 15%
random % weekday (Intercept) 0% 0% 0%
random % year (Intercept) 1% 0% 2%
random % Residual 55% 36% 73%

Continuous variables were z-transformed.

Exploration of Google Mobility

g_ <- fread(here::here("Data/google_mobility.txt")) # fwrite(d, here::here('Data/data.txt'), sep ='\t')
g_[, Year := as.integer(substring(date, nchar(date) - 3, nchar(date)))]
g_[nchar(date) == 9, date := paste0("0", date)]
g_[, date_ := as.Date(date, format = "%d.%m.%Y")]
g_[, Day := yday(date_)]
setnames(g_, old = "country_region", new = "Country")
g_[, weekday := weekdays(date_)]

g0 = ggplot(g_, aes(x = parks_percent_change_from_baseline, fill = factor(Year))) +
  geom_histogram(position = "dodge") +
  # scale_y_continuous(trans = 'log')+
  scale_fill_manual(values = c("orange", "skyblue", "black"), guide = 'none') +
  geom_vline(xintercept = 0, lty = 3, col = "red") +
  labs(subtitle = "Distribution") +
  xlab( 'Google Mobility\n[% change in human presence]') +
  ylab( 'Count') +
  facet_wrap(~Country, nrow = 5)

g1 = ggplot(g_, aes(x = Day, y = parks_percent_change_from_baseline, col = factor(Year))) +
  geom_line() +
  facet_wrap(~Country, nrow = 5) +
  labs(subtitle = "Raw data", xlab = 'Day\n ', y = 'Google Mobility\n[% change in human presence]') +
  # scale_y_continuous(trans = 'log')+
  coord_cartesian(ylim = c(-100, 300))+
  scale_color_manual(values = c("orange", "skyblue", "black"), guide = 'none')

g2 = ggplot(g_, aes(x = Day, y = parks_percent_change_from_baseline, col = factor(Year))) +
     stat_smooth() +
     facet_wrap(~Country, nrow = 5) +
     labs(subtitle = "Loess", xlab = 'Day\n ') +
     # scale_y_continuous(trans = 'log')+
     coord_cartesian(ylim = c(-100, 300))+
     scale_color_manual(values = c("orange", "skyblue", "black"), name = 'Year')+
      theme(
         axis.text.y = element_blank(),
         axis.title.y = element_blank()
       )
  ggarrange(
    g0, g1, g2,
    ncol = 3, widths = c(0.9, 1, 1.1)
  )

 ggsave(here::here("Outputs/Fig_ZZ.png"), width = 8*2.54, height = 5*2.54, unit = "cm", dpi = 600)

Fig. ZZ | Changes in human presence (Google Mobility) in parks across year and between years. Left plots represent the raw data, right plots LOESS smoothed curves. Google Mobility is absent for years before COVID-19. Nevertheless, 2022 was a year without shutdowns in the studied countries. Assuming that human activities levels might have been similar to pre-COVID-19 years (which may not be the case), human acctivity in the shutdown years (2020 and 2021) decreased. However, such decrease might be irrelevant for birds as the day-to-day variatioin in human presence seems larger than the general decrease in activity. For weekday specific plot see Fig. S_ZZ

g[, weekday := factor(weekday, levels = (c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday")))]

ggplot(g, aes(x = Day, y = parks_percent_change_from_baseline, col = factor(Year))) +
    geom_line() +
    facet_grid(rows = vars(Country), cols = vars(weekday)) +
    # scale_y_continuous(trans = 'log')+
    scale_color_manual(values = c("orange", "skyblue", "black"))

Fig. S_ZZ | Changes in human presence (Google Mobility) in parks across weekdays and years. Depicted are raw data.

Google Mobility vs Stringency

# Predictions 
l = list()
 sc = s[Country == "Czechia"]
 cz <- lmer(parks_percent_change_from_baseline ~
    StringencyIndex + 
   (scale(StringencyIndex)|weekday),
 # (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
 data = sc, REML = FALSE
 )
 bsim <- sim(cz, n.sim = nsim)
 v <- apply(bsim@fixef, 2, quantile, prob = c(0.5))
 newD <- data.frame(StringencyIndex = seq(min(sc$StringencyIndex), max(sc$StringencyIndex), length.out = 100)) # values to predict for
 X <- model.matrix(~StringencyIndex, data = newD) # exactly the model which was used has to be specified here
 newD$pred <- (X %*% v)
 predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
 for (j in 1:nsim) predmatrix[, j] <- (X %*% bsim@fixef[j, ])
 newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
 newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
 newD$pred <- apply(predmatrix, 1, quantile, prob = 0.5)
 newD$Country = 'Czechia'
 l[[1]] = newD
 
 s[, year_weekday :=paste(Year, weekday)]
 sf = s[Country == "Finland"]
 fi <- lmer(parks_percent_change_from_baseline ~
    Year+
    StringencyIndex + 
   (scale(StringencyIndex)|year_weekday),
 # (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
 data = sf, REML = FALSE
 )
 bsim <- sim(fi, n.sim = nsim)
 v <- apply(bsim@fixef, 2, quantile, prob = c(0.5))
 newD <- data.frame(Year = mean(sf$Year), StringencyIndex = seq(min(sf$StringencyIndex), max(sf$StringencyIndex), length.out = 100)) # values to predict for
 X <- model.matrix(~Year + StringencyIndex, data = newD) # exactly the model which was used has to be specified here
 newD$pred <- (X %*% v)
 predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
 for (j in 1:nsim) predmatrix[, j] <- (X %*% bsim@fixef[j, ])
 newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
 newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
 newD$pred <- apply(predmatrix, 1, quantile, prob = 0.5)
 newD$Country = 'Finland'
 newD$Year = NULL
 l[[2]] = newD

 sh <- s[Country == "Hungary"]
 hu <- lmer(parks_percent_change_from_baseline ~
     Year +
     StringencyIndex +
     (scale(StringencyIndex) | year_weekday),
 # (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
 data = sh, REML = FALSE
 )

 bsim <- sim(hu, n.sim = nsim)
 v <- apply(bsim@fixef, 2, quantile, prob = c(0.5))
 newD <- data.frame(Year = mean(sh$Year), StringencyIndex = seq(min(sh$StringencyIndex), max(sh$StringencyIndex), length.out = 100)) # values to predict for
 X <- model.matrix(~ Year + StringencyIndex, data = newD) # exactly the model which was used has to be specified here
 newD$pred <- (X %*% v)
 predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
 for (j in 1:nsim) predmatrix[, j] <- (X %*% bsim@fixef[j, ])
 newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
 newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
 newD$pred <- apply(predmatrix, 1, quantile, prob = 0.5)
 newD$Country <- "Hungary"
 newD$Year <- NULL
 l[[3]] <- newD
 
 sp <- s[Country == "Poland"]
 pl <- lmer(parks_percent_change_from_baseline ~
     Year +
     StringencyIndex +
     (scale(StringencyIndex) | year_weekday),
 # (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
 data = sp, REML = FALSE
 )
 bsim <- sim(pl, n.sim = nsim)
 v <- apply(bsim@fixef, 2, quantile, prob = c(0.5))
 newD <- data.frame(Year = mean(sp$Year), StringencyIndex = seq(min(sp$StringencyIndex), max(sp$StringencyIndex), length.out = 100)) # values to predict for
 X <- model.matrix(~ Year + StringencyIndex, data = newD) # exactly the model which was used has to be specified here
 newD$pred <- (X %*% v)
 predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
 for (j in 1:nsim) predmatrix[, j] <- (X %*% bsim@fixef[j, ])
 newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
 newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
 newD$pred <- apply(predmatrix, 1, quantile, prob = 0.5)
 newD$Country <- "Poland"
 newD$Year <- NULL
 l[[4]] <- newD

 sa <- s[Country == "Australia"]
 au <- lmer(parks_percent_change_from_baseline ~
     Year +
     StringencyIndex +
     (scale(StringencyIndex) | year_weekday),
 # (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
 data = sa, REML = FALSE
 )

 bsim <- sim(au, n.sim = nsim)
 v <- apply(bsim@fixef, 2, quantile, prob = c(0.5))
 newD <- data.frame(Year = mean(sa$Year), StringencyIndex = seq(min(sa$StringencyIndex), max(sa$StringencyIndex), length.out = 100)) # values to predict for
 X <- model.matrix(~ Year + StringencyIndex, data = newD) # exactly the model which was used has to be specified here
 newD$pred <- (X %*% v)
 predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
 for (j in 1:nsim) predmatrix[, j] <- (X %*% bsim@fixef[j, ])
 newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
 newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
 newD$pred <- apply(predmatrix, 1, quantile, prob = 0.5)
 newD$Country <- "Australia"
 newD$Year <- NULL
 l[[5]] <- newD

# Figure G_S
g_s = data.table(do.call(rbind,l))
g_s[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia")))]

col3_ = c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1]
col3__ = col3_[3:7]
#p = 
ggplot(g_s, aes(x = StringencyIndex, y = pred, col = Country)) +
  geom_ribbon(aes(ymin = lwr, ymax = upr, fill = Country, color = NULL), alpha = .15) +
  geom_jitter(aes(y = parks_percent_change_from_baseline, fill = Country), data = s, pch = 21, col = 'grey20', width = 0.7, height = 3, alpha = 0.5) +
  geom_line(lwd = 1) +
  labs(subtitle = "Mixed model per country predicitons", y = "Google Mobiligy\n[% change from baseline]", x = "Stringency Index") +
   # scale_color_locuszoom()+
   # scale_fill_locuszoom(guide = "none")
  scale_x_continuous(breaks = round(seq(25, 75, by = 25), 1)) +
  scale_y_continuous(breaks = round(seq(-100, 200, by = 50), 1)) +
  #scale_y_continuous(breaks = round(seq(-100, 175, by = 25), 1)) +
  scale_colour_manual(values = col3__, guide = guide_legend(reverse = TRUE, override.aes = list(size = 0)), 
            labels = paste("<span style='color:",
                                   col3__,
                                   "'>",
                                   levels(sp$Country),
                                   "</span>")
            ) +
  scale_fill_manual(values = col3__, guide = "none") +
  theme_bw() +
  theme(
    legend.text = element_markdown(size = 6),
    #legend.position = "right",
    legend.title = element_blank(),
    # legend.spacing.y = unit(0.1, 'cm'),
    legend.key.height = unit(0.5, "line"),
    legend.key.size = unit(0, "line"),
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin=margin(-10,1,-10,-10),
    # legend.position=c(0.5,1.6),
    plot.title = element_text(color = "grey", size = 7),
    plot.subtitle = element_text(color = "grey60", size = 6),
    plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 0.5, unit = "pt"),
    panel.grid = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(colour = ax_lines, size = 0.25),
    axis.ticks = element_line(colour = ax_lines, size = 0.25),
    # axis.text.x = element_text()
    axis.ticks.length = unit(1, "pt"),
    axis.text = element_text(, size = 6),
    axis.title = element_text(size = 7)
  )

ggsave(here::here("Outputs/Fig_G-S_rev_width_CustomLocusZoom_year-weekday.png"), width = 7, height = 6, unit = "cm", dpi = 600)

Fig. 2 | Relationship betweeen Google Mobility and Stringency Index. Lines with shaded areas represent predicted relationships from country-specific mixed effect models controlled for the year and non-independence of data points by including weekday within the year as random intercept and Stringency Index as a random slope. Dots represent raw data, jittered to increase visibility. Colors indicate country. The predictions reveal generally negative and week relationship between Gooble Mobility (human activity) in parks and seveareness of governmental restrictions.

Table S_G_S | Google Mobility in relations to stringency index

ll = list()
s[, year_weekday := paste(Year, weekday)]

sf = s[Country == "Finland"]
fi <- lmer(scale(parks_percent_change_from_baseline) ~
  scale(Year) +
  scale(StringencyIndex) +
  (scale(StringencyIndex) | year_weekday),
data = sf, REML = FALSE
)
ll[[1]] = m_out(name = "Table S_g-s - FI", dep = "Google Mobility", model = fi, nsim = 5000)

sp <- s[Country == "Poland"]
pl <- lmer(scale(parks_percent_change_from_baseline) ~
  scale(Year) +
  scale(StringencyIndex) +
  (scale(StringencyIndex) | year_weekday),
data = sp, REML = FALSE
)
ll[[2]] = m_out(name = "Table S_g-s - PL", dep = "Google Mobility", model = pl, nsim = 5000)

sc = s[Country == "Czechia"]
cz <- lmer(scale(parks_percent_change_from_baseline) ~
  scale(StringencyIndex) +
  (scale(StringencyIndex) | weekday),
data = sc, REML = FALSE
)
ll[[3]] = m_out(name = "Table S_g-s - CZ", dep = "Google Mobility", model = cz, nsim = 5000)

sh <- s[Country == "Hungary"]
hu <- lmer(scale(parks_percent_change_from_baseline) ~
  scale(Year) +
  scale(StringencyIndex) +
  (scale(StringencyIndex) | year_weekday),
data = sh, REML = FALSE
)
ll[[4]] = m_out(name = "Table S_g-s - HU", dep = "Google Mobility", model = hu, nsim = 5000)

sa <- s[Country == "Australia"]
au <- lmer(scale(parks_percent_change_from_baseline) ~
  scale(Year) +
  scale(StringencyIndex) +
  (scale(StringencyIndex) | year_weekday),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = sa, REML = FALSE
)
ll[[5]] = m_out(name = "Table S_g-s - AU", dep = "Google Mobility", model = au, nsim = 5000)

out_g_s = data.table(do.call(rbind, ll))
out_g_s[is.na(out_g_s)] <- ""
out_g_s$R2_mar = out_g_s$R2_con = NULL
out_g_s[, effect := gsub("scale\\(Year\\)", "year", effect)]
out_g_s[, effect := gsub("scale\\(StringencyIndex\\)", "stringency index", effect)]
out_g_s[, effect := gsub("year_weekday", "weekday within year", effect)]
fwrite(file = here::here("Outputs/Table_S_G-S.csv"), out_g_s)

out_g_s %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
model response error_structure N type effect estimate_r lwr_r upr_r
Table S_g-s - FI Google Mobility Gaussian 322 fixed (Intercept) -0.268 -1.064 0.503
fixed year 0.939 0.459 1.424
fixed stringency index -0.259 -1.004 0.494
random % weekday within year (Intercept) 22% 47% 54%
random % weekday within year stringency index 22% 47% 54%
random % Residual 55% -7% 5%
Table S_g-s - PL Google Mobility Gaussian 329 fixed (Intercept) 0.013 -0.605 0.637
fixed year 0.78 0.011 1.535
fixed stringency index 0.025 -0.362 0.412
random % weekday within year (Intercept) 47% 44% 49%
random % weekday within year stringency index 47% 44% 49%
random % Residual 6% 2% 13%
Table S_g-s - CZ Google Mobility Gaussian 1054 fixed (Intercept) 0.006 -0.526 0.536
fixed stringency index -0.219 -0.635 0.189
random % weekday (Intercept) 23% -184% 40%
random % weekday stringency index 23% -184% 40%
random % Residual 54% 21% 469%
Table S_g-s - HU Google Mobility Gaussian 874 fixed (Intercept) 0.255 -0.272 0.737
fixed year 0.007 -0.457 0.466
fixed stringency index -1.264 -1.946 -0.566
random % weekday within year (Intercept) 48% 2% 49%
random % weekday within year stringency index 48% 2% 49%
random % Residual 4% 2% 95%
Table S_g-s - AU Google Mobility Gaussian 1065 fixed (Intercept) 0.528 -0.01 1.074
fixed year 0.441 0.245 0.638
fixed stringency index -0.54 -0.855 -0.22
random % weekday within year (Intercept) -31% 44% 65%
random % weekday within year stringency index -31% 44% 65%
random % Residual 162% -30% 13%

Continuous variables were z-transformed.

FID ~ Stringency

# predictions for fig
 # full
 mss <- lmer(scale(log(FID)) ~
     scale(Year) +
     scale(log(SD)) +
     scale(log(FlockSize)) +
     scale(log(BodyMass)) +
     scale(sin(rad)) + scale(cos(rad)) +
     # scale(Day)+
     scale(Temp) +
     scale(StringencyIndex) +
     (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1 | IDLocality) + (1 | sp_loc),
 data = s, REML = FALSE,
 control = lmerControl(
     optimizer = "optimx", optCtrl = list(method = "nlminb")
 )
 )
 est_mss <- est_out(mss, "ALL: (1|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc)")
 est_mss[, control_for_starting_distance := "yes"]

 msx <- lmer(scale(log(FID)) ~
    scale(Year) +
    #scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1| IDLocality) + (1 | sp_loc),
 data = s, REML = FALSE,
 control = lmerControl(
    optimizer = "optimx", optCtrl = list(method = "nlminb")
)
)
est_msx <- est_out(msx, "ALL: (1|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc)")
est_msx[, control_for_starting_distance := "no"]


# CZ - singular fits only due to genera estimated as zero (removing it changes no results)
  css <- lmer(scale(log(FID)) ~
    #scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
    #(1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
    data = s[Country == "Czechia"], REML = FALSE
    )
  est_css <- est_out(css, "Czechia: (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
  est_css[, control_for_starting_distance := "yes"]
  
  csx <- lmer(scale(log(FID)) ~
    #scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
    #(1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
    data = s[Country == "Czechia"], REML = FALSE
    )
  est_csx <- est_out(csx, "Czechia: (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
  est_csx[, control_for_starting_distance := "no"]

# FI
fss <- lmer(scale(log(FID)) ~
    scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
data = s[Country == "Finland"], REML = FALSE
)
est_fss <- est_out(fss, "Finland: (1|genus)+(1|Species)+(1|sp_day_year)+(scale(StringencyIndex)|IDLocality)+(1|sp_loc)")
est_fss[, control_for_starting_distance := "yes"]

fsx <- lmer(scale(log(FID)) ~
    scale(Year) +
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
data = s[Country == "Finland"], REML = FALSE
)
est_fsx <- est_out(fsx, "Finland: (1|genus)+(1|Species)+(1|sp_day_year)+(scale(StringencyIndex)|IDLocality)+(1|sp_loc)")
est_fsx[, control_for_starting_distance := "no"]

# HU - singular fits only due to sp_loc estimated as zero (removing it changes no results)
hss <- lmer(scale(log(FID)) ~
    scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
     (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Hungary"], REML = FALSE
)
est_hss <- est_out(hss, "Hungary: (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_hss[, control_for_starting_distance := "yes"]

hsx <- lmer(scale(log(FID)) ~
    scale(Year) +
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Hungary"], REML = FALSE
)
est_hsx <- est_out(hsx, "Hungary: (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_hsx[, control_for_starting_distance := "no"]

# AU - singular fits only due to Year and random slope estimated as zero (removing those changes no results)
ass <- lmer(scale(log(FID)) ~
    scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
      (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Australia"], REML = FALSE
)
est_ass <- est_out(ass, "Australia: (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_ass[, control_for_starting_distance := "yes"]

asx <- lmer(scale(log(FID)) ~
    scale(Year) +
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Australia"], , REML = FALSE,
        control = lmerControl( 
            optimizer ='optimx', optCtrl=list(method='nlminb')) 
)
est_asx <- est_out(asx, "Australia: (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_asx[, control_for_starting_distance := "no"]

# PL 
pss <- lmer(scale(log(FID)) ~
    scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1 | genus) + (1 | Species) + (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Poland"], REML = FALSE
)
est_pss <- est_out(pss, "Poland: (1|genus)+(1|Species)+(1|sp_day_year)")
est_pss[, control_for_starting_distance := "yes"]

psx <- lmer(scale(log(FID)) ~
    scale(Year) +
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1 | genus) + (1 | Species) + (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Poland"], REML = FALSE
)
est_psx <- est_out(psx, "Poland: (1|genus)+(1|Species)+(1|sp_day_year)")
est_psx[, control_for_starting_distance := "no"]

  # combine
    est_mss[, Country := 'All\n(mixed model)']
    est_msx[, Country := "All\n(mixed model)"]
    est_ass[, Country := "Australia"]
    est_asx[, Country := "Australia"]
    est_css[, Country := "Czechia"]
    est_csx[, Country := "Czechia"]
    est_hss[, Country := "Hungary"]
    est_hsx[, Country := "Hungary"]
    est_pss[, Country := "Poland"]
    est_psx[, Country := "Poland"]
    est_fss[, Country := "Finland"]
    est_fsx[, Country := "Finland"]

    os = rbind(est_mss, est_msx, 
            est_ass, est_asx, 
            est_css, est_csx, 
            est_hss, est_hsx,
            est_pss, est_psx, 
            est_fss, est_fsx)
    save(os, file = here::here('Data/dat_est_Stringency_rev.Rdata'))

# estimates for table
  mss_out <- m_out(name = "Table S_FID_s - full a", dep = "Escape distance", model = mss, nsim = 5000)
  msx_out <- m_out(name = "Table S_FID_s - full b", dep = "Escape distance", model = msx, nsim = 5000)
  css_out <- m_out(name = "Table v - CZ a", dep = "Escape distance", model = css, nsim = 5000)
  csx_out <- m_out(name = "Table S_FID_s - CZ b", dep = "Escape distance", model = csx, nsim = 5000)
  fss_out <- m_out(name = "Table S_FID_s - FI a", dep = "Escape distance", model = fss, nsim = 5000)
  fsx_out <- m_out(name = "Table S_FID_s - FI b", dep = "Escape distance", model = fsx, nsim = 5000)
  hss_out <- m_out(name = "Table S_FID_s - HU a", dep = "Escape distance", model = hss, nsim = 5000)
  hsx_out <- m_out(name = "Table S_FID_s - HU b", dep = "Escape distance", model = hsx, nsim = 5000)
  ass_out <- m_out(name = "Table S_FID_s - AU a", dep = "Escape distance", model = ass, nsim = 5000)
  asx_out <- m_out(name = "Table S_FID_s - AU b", dep = "Escape distancey", model = asx, nsim = 5000)
  pss_out <- m_out(name = "Table S_FID_s - PL a", dep = "Escape distance", model = pss, nsim = 5000)
  psx_out <- m_out(name = "Table S_FID_s - PL b", dep = "Escape distancey", model = psx, nsim = 5000)

  out_FID_s <- rbind(mss_out, msx_out, fss_out, fsx_out, pss_out, psx_out, css_out, csx_out, hss_out, hsx_out, ass_out, asx_out, fill = TRUE)
  out_FID_s[is.na(out_FID_s)] <- ""
  out_FID_s$R2_mar = out_FID_s$R2_con = NULL
  out_FID_s[, effect := gsub("scale\\(StringencyIndex\\)", "Stringency index", effect)]
  out_FID_s[, effect := gsub("scale\\(Year\\)", "year", effect)]
  out_FID_s[, effect := gsub("scale\\(log\\(SD\\)", "starting distance (ln)", effect)]
  out_FID_s[, effect := gsub("scale\\(Temp\\)", "temperaturre", effect)]
  out_FID_s[, effect := gsub("scale\\(log\\(FlockSize\\)\\)", "flock size", effect)]
  out_FID_s[, effect := gsub("scale\\(log\\(BodyMass\\)\\)", "body mass", effect)]
  out_FID_s[, effect := gsub("scale\\(sin\\(rad\\)\\)", "sine of radians", effect)]
  out_FID_s[, effect := gsub("scale\\(cos\\(rad\\)\\)", "cosine of radians", effect)]
  out_FID_s[, effect := gsub("Species", "species", effect)]
  out_FID_s[, effect := gsub("sp_day_year", "species within day & year", effect)]
  out_FID_s[, effect := gsub("IDLocality", "site", effect)]
  out_FID_s[, effect := gsub("sp_loc", "species within site", effect)]
  fwrite(file = here::here("Outputs/Table_S_FID_s_c.csv"), out_FID_s)
load(here::here("Data/dat_est_Stringency_rev.Rdata"))
os[predictor %in% c("scale(StringencyIndex)"), predictor := "Stringency Index"]
oso <- os[predictor %in% c("Stringency Index")]
oso[, N:=as.numeric(sub('.*N = ', '', model))]
# add meta-analytical mean
  oso_s = oso[control_for_starting_distance == 'yes']
  met = summary(meta.summaries(d = oso_s$estimate, se = oso_s$sd, method = "fixed", weights = oso_s$N))$summci
  oso_met = data.table(predictor = "Period", estimate = met[2], lwr = met[1], upr = met[3], sd = NA, model = NA, control_for_starting_distance = "yes", Country = "Combined\n(metanalytical)", N = NA)

  oso_sx = oso[control_for_starting_distance == "no"]
  metx = summary(meta.summaries(d = oso_sx$estimate, se = oso_sx$sd, method = "fixed", weights = oso_sx$N))$summci
  oso_metx = data.table(predictor = "Period", estimate = metx[2], lwr = metx[1], upr = metx[3], sd = NA, model = NA, control_for_starting_distance = "no", Country = "Combined\n(metanalytical)", N = NA)
  
  oso = rbind(oso, oso_met, oso_metx)
    
oso[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia", "Combined\n(metanalytical)", "All\n(mixed model)")))]

width_ <- .5 # spacing between error bars

#col_ <- c(brewer.pal(n = 12, name = "Paired"), "grey30", "grey80")
#Tol_bright <- c("#EE6677", "#228833", "#4477AA", "#CCBB44", "#66CCEE", "#AA3377", "#BBBBBB")
#Tol_muted <- c("#88CCEE", "#44AA99", "#117733", "#332288", "#DDCC77", "#999933", "#CC6677", "#882255", "#AA4499", "#DDDDDD")
#Tol_light <- c("#BBCC33", "#AAAA00", "#77AADD", "#EE8866", "#EEDD88", "#FFAABB", "#99DDFF", "#44BB99", "#DDDDDD")

# From Color Universal Design (CUD): https://jfly.uni-koeln.de/color/
#Okabe_Ito <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#000000")
#col_ = Okabe_Ito[7:1]
# JAMA and LocusZoom modified order
#col_ =  c("#374E55FF", "#374E55FF", "#DF8F44FF", "#79AF97FF", "#00A1D5FF", "#B24745FF",  "#80796BFF") #"#6A6599FF",
#col_ <- c("#357EBDFF", "#9632B8FF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#D43F3AFF", "#D43F3AFF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
col_ = c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
#show_col(col_)

#p = 
ggplot(oso, aes(x = estimate, y = Country, col = Country, shape = control_for_starting_distance)) +
    geom_vline(xintercept = 0, color = "grey", linetype = "dotted") +
    geom_errorbarh(aes(xmin = lwr, xmax = upr), height = 0, position = ggstance::position_dodgev(width_)) +
    # geom_point(position = ggstance::position_dodgev(.6)) +
    geom_point(position = position_dodge(width = width_), bg = "white", size = 1.1) +
    # scale_color_viridis(discrete=TRUE, begin=0, end = 0.5)  +
    # scale_fill_viridis(discrete=TRUE, begin=0, end = 0.5) +
    # geom_text( aes(x = n_pos,label = N), vjust = 0, size = 1.75, position = ggstance::position_dodgev(width_))+ # 3 positions for 3 bars
    # annotate("text", x=log10(3), y=85, label= "Used", col = "grey30", size = 2.5)+

    scale_shape_manual(name = "Controlled for\nstarting distance", guide = guide_legend(reverse = TRUE), values = c(21, 19)) +
    #scale_color_jama(guide = "none")+ #, palette = 'light'
    scale_color_manual(guide = "none", values = col_) + #guide_legend(reverse = TRUE)
    scale_x_continuous(breaks = round(seq(-0.3, 0.4, by = 0.1), 1)) +
    ylab("") +
    xlab("Standardized effect size of Stringency Index\n[on flight initiation distance]") +
    # coord_cartesian(xlim = c(-.15, .15)) +
    # scale_x_continuous(breaks = round(seq(-.15, .15, by = 0.05),2)) +
    theme_bw() +
    theme(
        legend.position = "right",
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 6),
        # legend.spacing.y = unit(0.1, 'cm'),
        legend.key.height = unit(0.5, "line"),
        legend.margin = margin(0, 0, 0, 0),
        # legend.position=c(0.5,1.6),
        plot.title = element_text(color = "grey", size = 7),
        plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 0.5, unit = "pt"),
        panel.grid = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = ax_lines, size = 0.25),
        axis.line.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
        # axis.text.x = element_text()
        axis.ticks.length = unit(1, "pt"),
        axis.text.x = element_text(, size = 6),
        axis.text.y = element_text(colour = "black", size = 7),
        axis.title = element_text(size = 7)
    )

ggsave(here::here("Outputs/Fig_Stringency_rev_width_CustomLocusZoom_v2.png"), width = 8, height = 6, unit = "cm", dpi = 600)

Fig. S_stringency | The effect size with 95%CIs of Stringency Index on flight initiation distance (ln-transformed). The general model on all data was controlled for year, starting distance (ln-transformed; filled circles) or not (empty circles), flock size (ln-transformed), temperature (also a proxy for a day within the breeding season: r Pearson = 0.48; Fig. S1) and time of day. To account for circular properties of time, time was transformed into radians (2 × time × π/24) and fitted as sine and cosine of radians (Bulla et al. 2016). All continuous variables were standardised by subtracting the mean and dividing by the standard deviation. The multicollinearity was small as correlations between predictors were weak (r<XX, Fig. S1). To account for the non-independence of data points (Schielzeth & Forstmeier 2009; Barr et al. 2013), we fitted random intercepts of genus, species, species at a given day and year, site, and species within a site and in the full model Strigency also random slope of Stringency within genus and country. In addition, we attempted to fit Stringency as random slope within genus for country-specific modl, but it explained little variance and estimates remained the same. Note that in the full model, fitting Stringency as random slope at other random intercepts produces similar results (see Fig Sxx). We used this approach with a full dataset with all observations (n = xx), as well as with conservative datasets, one with at least five observations per species , the other with at least 10 observations per species. Albeit random structure of the full model accounts for the potential country specific effect (country accounted for xx% of variance, Table xx), depicted are also estimates from country specific models, with same predictors and random structure as the full model (excluding the random slopes, and random intercept of country and in case of Poland also site, as specific sites were not noted in Poland). Using meta.summaries function from rmeta R-package we used the country estimates, their standard deviation, and sample size per country to estimate the meta-analytical mean, which reflects the results based on the full mixed effect model. Importantly, the models test for within year changes in esccape distancces due to human changes presence, i.e. week to week changes in Stringency Index. In other words, the models investigate whether birds flexibly adjust their fear response, a test different to the one intended in this study. Albeit the response of humans to governmental restrictions were likely country specific and presence of humans in parks might have increased in some countries, decreased in others or did not change, the escape distances do not strongly reflect such changes. In other words, the birds are either inflexible or the change in human behavior due to shutdowns was not strong enought, which might have been the case - see Fig. ZZ - decide what to show in this figure (I think 2020-2022 changes in Google Mobility for each country including the fits and refering to supplementary figure). Note however that human presence consistently declined with increased restrictions across countries, albeit the relationships are weak (Fig. 2).

Table FID_s | Escape distance in relations to stringency index

out_FID_s %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
model response error_structure N type effect estimate_r lwr_r upr_r
Table S_FID_s - full a Escape distance Gaussian 3676 fixed (Intercept) -0.076 -0.334 0.187
fixed year 0.02 -0.027 0.066
fixed starting distance (ln)) 0.491 0.46 0.523
fixed flock size -0.003 -0.03 0.024
fixed body mass 0.023 -0.047 0.092
fixed sine of radians -0.02 -0.061 0.023
fixed cosine of radians -0.005 -0.04 0.03
fixed temperaturre -0.007 -0.047 0.035
fixed Stringency index 0.015 -0.157 0.183
random % species within day & year (Intercept) 5% 4% 7%
random % species within site (Intercept) 5% 4% 6%
random % species (Intercept) 8% 6% 8%
random % site (Intercept) 8% 7% 9%
random % genus (Intercept) 4% 4% 4%
random % Country (Intercept) 1% -14% 14%
random % Country Stringency index 1% -14% 14%
random % Residual 68% 47% 93%
Table S_FID_s - full b Escape distance Gaussian 3676 fixed (Intercept) -0.113 -0.458 0.241
fixed year 0.085 0.028 0.141
fixed flock size 0.01 -0.019 0.04
fixed body mass 0.136 0.049 0.223
fixed sine of radians 0.045 -0.002 0.093
fixed cosine of radians 0.048 0.008 0.088
fixed temperaturre 0.008 -0.039 0.055
fixed Stringency index 0.01 -0.216 0.234
random % species within day & year (Intercept) 6% 4% 10%
random % species within site (Intercept) 5% 4% 9%
random % species (Intercept) 10% 7% 14%
random % site (Intercept) 11% 8% 15%
random % genus (Intercept) 6% 5% 8%
random % Country (Intercept) 1% -26% 18%
random % Country Stringency index 1% -26% 18%
random % Residual 58% 36% 97%
Table S_FID_s - FI a Escape distance Gaussian 354 fixed (Intercept) 0.033 -0.179 0.243
fixed year -0.105 -0.337 0.121
fixed starting distance (ln)) 0.187 0.086 0.291
fixed flock size -0.111 -0.21 -0.008
fixed body mass 0.177 0.012 0.34
fixed sine of radians 0.06 -0.062 0.184
fixed cosine of radians 0.012 -0.109 0.135
fixed temperaturre -0.084 -0.197 0.028
fixed Stringency index 0.156 -0.074 0.381
random % species within site (Intercept) 12% 12% 13%
random % species within day & year (Intercept) 3% 3% 3%
random % site (Intercept) 12% 10% 15%
random % species (Intercept) 0% 0% 1%
random % genus (Intercept) 3% 2% 5%
random % Residual 68% 64% 73%
Table S_FID_s - FI b Escape distance Gaussian 354 fixed (Intercept) -0.012 -0.231 0.209
fixed year -0.067 -0.309 0.167
fixed flock size -0.107 -0.206 -0.003
fixed body mass 0.249 0.079 0.414
fixed sine of radians 0.054 -0.074 0.179
fixed cosine of radians 0.005 -0.117 0.124
fixed temperaturre -0.086 -0.2 0.025
fixed Stringency index 0.137 -0.092 0.374
random % species within site (Intercept) 14% 14% 15%
random % species within day & year (Intercept) 4% 3% 4%
random % site (Intercept) 14% 11% 16%
random % species (Intercept) 3% 2% 5%
random % genus (Intercept) 1% 1% 2%
random % Residual 64% 59% 69%
Table S_FID_s - PL a Escape distance Gaussian 329 fixed (Intercept) 0.002 -0.128 0.136
fixed year -0.161 -0.394 0.073
fixed starting distance (ln)) 0.532 0.447 0.617
fixed flock size -0.024 -0.104 0.06
fixed body mass -0.098 -0.205 0.009
fixed sine of radians 0.102 -0.16 0.36
fixed cosine of radians -0.157 -0.416 0.112
fixed temperaturre -0.171 -0.276 -0.065
fixed Stringency index 0.027 -0.207 0.262
random % species within day & year (Intercept) 13% 12% 14%
random % species (Intercept) 8% 6% 10%
random % genus (Intercept) 0% 0% 0%
random % Residual 79% 77% 82%
Table S_FID_s - PL b Escape distancey Gaussian 329 fixed (Intercept) 0.001 -0.199 0.203
fixed year -0.363 -0.634 -0.101
fixed flock size 0.06 -0.037 0.16
fixed body mass -0.028 -0.181 0.121
fixed sine of radians -0.006 -0.304 0.305
fixed cosine of radians -0.043 -0.344 0.276
fixed temperaturre -0.162 -0.288 -0.037
fixed Stringency index -0.027 -0.298 0.24
random % species within day & year (Intercept) 7% 7% 7%
random % species (Intercept) 22% 17% 26%
random % genus (Intercept) 0% 0% 0%
random % Residual 72% 67% 77%
Table v - CZ a Escape distance Gaussian 1054 fixed (Intercept) 0.223 0.025 0.437
fixed starting distance (ln)) 0.354 0.295 0.413
fixed flock size 0.004 -0.045 0.056
fixed body mass 0.181 0.026 0.341
fixed sine of radians 0.062 -0.013 0.136
fixed cosine of radians 0.007 -0.063 0.077
fixed temperaturre 0.011 -0.076 0.098
fixed Stringency index -0.062 -0.181 0.058
random % species within site (Intercept) 4% 4% 4%
random % species within day & year (Intercept) 2% 2% 2%
random % species (Intercept) 15% 12% 17%
random % genus (Intercept) 14% 10% 18%
random % site (Intercept) 5% 3% 6%
random % Residual 61% 54% 68%
Table S_FID_s - CZ b Escape distance Gaussian 1054 fixed (Intercept) 0.229 0.02 0.433
fixed starting distance (ln)) 0.355 0.297 0.413
fixed flock size 0.004 -0.045 0.056
fixed body mass 0.182 0.028 0.342
fixed sine of radians 0.061 -0.012 0.137
fixed cosine of radians 0.006 -0.064 0.078
fixed temperaturre 0.01 -0.08 0.098
fixed Stringency index -0.061 -0.179 0.056
random % species within site (Intercept) 4% 4% 4%
random % species within day & year (Intercept) 2% 2% 2%
random % species (Intercept) 15% 12% 17%
random % genus (Intercept) 14% 10% 18%
random % site (Intercept) 5% 3% 6%
random % Residual 61% 54% 68%
Table S_FID_s - HU a Escape distance Gaussian 874 fixed (Intercept) 0.069 -0.212 0.342
fixed year 0.193 0.089 0.301
fixed starting distance (ln)) 0.507 0.447 0.569
fixed flock size 0.002 -0.052 0.054
fixed body mass -0.001 -0.14 0.141
fixed sine of radians -0.124 -0.19 -0.057
fixed cosine of radians 0.033 -0.039 0.105
fixed temperaturre -0.008 -0.103 0.09
fixed Stringency index 0.045 -0.096 0.183
random % species within day & year (Intercept) 3% 3% 3%
random % species within site (Intercept) 0% 0% 0%
random % species (Intercept) 1% 1% 1%
random % genus (Intercept) 11% 7% 15%
random % site (Intercept) 10% 7% 13%
random % Residual 74% 67% 82%
Table S_FID_s - HU b Escape distance Gaussian 874 fixed (Intercept) 0.055 -0.402 0.521
fixed year 0.236 0.107 0.365
fixed flock size 0.004 -0.056 0.063
fixed body mass 0.295 0.066 0.521
fixed sine of radians -0.14 -0.221 -0.064
fixed cosine of radians 0.034 -0.051 0.115
fixed temperaturre -0.031 -0.147 0.08
fixed Stringency index 0.003 -0.159 0.161
random % species within day & year (Intercept) 5% 5% 6%
random % species within site (Intercept) 0% 0% 0%
random % species (Intercept) 1% 1% 1%
random % genus (Intercept) 25% 18% 31%
random % site (Intercept) 18% 16% 19%
random % Residual 51% 43% 60%
Table S_FID_s - AU a Escape distance Gaussian 1065 fixed (Intercept) -0.049 -0.189 0.09
fixed year 0.036 -0.03 0.1
fixed starting distance (ln)) 0.496 0.444 0.549
fixed flock size 0.024 -0.024 0.072
fixed body mass -0.052 -0.146 0.044
fixed sine of radians -0.037 -0.109 0.035
fixed cosine of radians -0.052 -0.119 0.016
fixed temperaturre 0.033 -0.028 0.094
fixed Stringency index 0.058 -0.006 0.123
random % species within day & year (Intercept) 6% 6% 6%
random % species within site (Intercept) 12% 12% 12%
random % species (Intercept) 11% 9% 13%
random % genus (Intercept) 2% 2% 3%
random % site (Intercept) 7% 5% 9%
random % Residual 61% 57% 66%
Table S_FID_s - AU b Escape distancey Gaussian 1065 fixed (Intercept) -0.101 -0.287 0.087
fixed year 0.104 0.029 0.176
fixed flock size 0.055 0.001 0.111
fixed body mass 0.062 -0.063 0.184
fixed sine of radians -0.011 -0.095 0.072
fixed cosine of radians -0.019 -0.095 0.059
fixed temperaturre 0.041 -0.029 0.11
fixed Stringency index 0.091 0.015 0.164
random % species within day & year (Intercept) 7% 7% 8%
random % species within site (Intercept) 11% 11% 11%
random % species (Intercept) 13% 11% 15%
random % genus (Intercept) 7% 6% 8%
random % site (Intercept) 8% 6% 10%
random % Residual 54% 49% 58%

Continuous variables were z-transformed.

FID ~ Google Mobiliity

# predictions for fig
 # full
 mgs <- lmer(scale(log(FID)) ~
     scale(Year) +
     scale(log(SD)) +
     scale(log(FlockSize)) +
     scale(log(BodyMass)) +
     scale(sin(rad)) + scale(cos(rad)) +
     # scale(Day)+
     scale(Temp) +
     scale(parks_percent_change_from_baseline) +
     (scale(parks_percent_change_from_baseline)| genus) + (1 | Species) + (1 | sp_day_year) + 
     (scale(parks_percent_change_from_baseline)| Country) + (1 | IDLocality) + (1 | sp_loc),
 data = ss, REML = FALSE,
 control = lmerControl(
     optimizer = "optimx", optCtrl = list(method = "nlminb")
 )
 )
 est_mgs <- est_out(mgs, "ALL: (scale(Google)|genus) + (1|Species)  + (1|sp_day_year) + (scale(Google)|Country) + (1|IDLocality) +(1|sp_loc)")
 est_mgs[, control_for_starting_distance := "yes"]

 mgx <- lmer(scale(log(FID)) ~
    scale(Year) +
    #scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    scale(parks_percent_change_from_baseline) +
     (scale(parks_percent_change_from_baseline)| genus) + (1 | Species) + (1 | sp_day_year) + 
     (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality) + (1 | sp_loc),
 data = ss, REML = FALSE,
 control = lmerControl(
    optimizer = "optimx", optCtrl = list(method = "nlminb")
)
)
est_mgx <- est_out(mgx, "ALL: (scale(Google)|genus) + (1|Species) + (1|sp_day_year) + (scale(Google)|Country) + (1|IDLocality) +(1|sp_loc)")
est_mgx[, control_for_starting_distance := "no"]


# CZ - singular fits only due to genera estimated as zero (removing it changes no results)
  cgs <- lmer(scale(log(FID)) ~
    #scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
    (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
    #(1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
    data = ss[Country == "Czechia"], REML = FALSE
    )
  est_cgs <- est_out(cgs, "Czechia: (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
  est_cgs[, control_for_starting_distance := "yes"]
  
  cgx <- lmer(scale(log(FID)) ~
    #scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
    (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
    #(1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
    data = ss[Country == "Czechia"], REML = FALSE
    )
  est_cgx <- est_out(cgx, "Czechia: (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
  est_cgx[, control_for_starting_distance := "no"]

# FI
fgs <- lmer(scale(log(FID)) ~
    scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
    (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
data = ss[Country == "Finland"], REML = FALSE
)
est_fgs <- est_out(fgs, "Finland: (1|genus)+(1|Species)+(1|sp_day_year)+(scale(StringencyIndex)|IDLocality)+(1|sp_loc)")
est_fgs[, control_for_starting_distance := "yes"]

fgx <- lmer(scale(log(FID)) ~
    scale(Year) +
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
    (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
data = ss[Country == "Finland"], REML = FALSE
)
est_fgx <- est_out(fgx, "Finland: (1|genus)+(1|Species)+(1|sp_day_year)+(scale(StringencyIndex)|IDLocality)+(1|sp_loc)")
est_fgx[, control_for_starting_distance := "no"]

# HU - singular fits only due to sp_loc estimated as zero (removing it changes no results)
hgs <- lmer(scale(log(FID)) ~
    scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
     (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Hungary"], REML = FALSE
)
est_hgs <- est_out(hgs, "Hungary: (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_hgs[, control_for_starting_distance := "yes"]

hgx <- lmer(scale(log(FID)) ~
    scale(Year) +
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
    (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Hungary"], REML = FALSE
)
est_hgx <- est_out(hgx, "Hungary: (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_hgx[, control_for_starting_distance := "no"]

# AU - singular fits only due to Year and random slope estimated as zero (removing those changes no results)
ags <- lmer(scale(log(FID)) ~
    scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
      (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Australia"], REML = FALSE
)
est_ags <- est_out(ags, "Australia: (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_ags[, control_for_starting_distance := "yes"]

agx <- lmer(scale(log(FID)) ~
    scale(Year) +
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
    (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Australia"], , REML = FALSE,
        control = lmerControl( 
            optimizer ='optimx', optCtrl=list(method='nlminb')) 
)
est_agx <- est_out(agx, "Australia: (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_agx[, control_for_starting_distance := "no"]

# PL 
pgs <- lmer(scale(log(FID)) ~
    scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
    (1 | genus) + (1 | Species) + (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Poland"], REML = FALSE
)
est_pgs <- est_out(pgs, "Poland: (1|genus)+(1|Species)+(1|sp_day_year)")
est_pgs[, control_for_starting_distance := "yes"]

pgx <- lmer(scale(log(FID)) ~
    scale(Year) +
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
    (1 | genus) + (1 | Species)+ (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Poland"], REML = FALSE
)
est_pgx <- est_out(pgx, "Poland: (1|genus)+(1|Species)+(1|sp_day_year)")
est_pgx[, control_for_starting_distance := "no"]

  # combine
    est_mgs[, Country := 'All\n(mixed model)']
    est_mgx[, Country := "All\n(mixed model)"]
    est_ags[, Country := "Australia"]
    est_agx[, Country := "Australia"]
    est_cgs[, Country := "Czechia"]
    est_cgx[, Country := "Czechia"]
    est_hgs[, Country := "Hungary"]
    est_hgx[, Country := "Hungary"]
    est_pgs[, Country := "Poland"]
    est_pgx[, Country := "Poland"]
    est_fgs[, Country := "Finland"]
    est_fgx[, Country := "Finland"]

    og = rbind(est_mgs, est_mgx, 
            est_ags, est_agx, 
            est_cgs, est_cgx, 
            est_hgs, est_hgx,
            est_pgs, est_pgx, 
            est_fgs, est_fgx)
    save(og, file = here::here('Data/dat_est_Google_rev.Rdata'))
  # estimatees for table
  mgs_out <- m_out(name = "Table S_FID_g - full a", dep = "Escape distance", model = mgs, nsim = 5000)
  mgx_out <- m_out(name = "Table S_FID_g - full b", dep = "Escape distance", model = mgx, nsim = 5000)
  cgs_out <- m_out(name = "Table S_FID_g - CZ a", dep = "Escape distance", model = cgs, nsim = 5000)
  cgx_out <- m_out(name = "Table S_FID_g - CZ b", dep = "Escape distance", model = cgx, nsim = 5000)
  fgs_out <- m_out(name = "Table S_FID_g - FI a", dep = "Escape distance", model = fgs, nsim = 5000)
  fgx_out <- m_out(name = "Table S_FID_g - FI b", dep = "Escape distance", model = fgx, nsim = 5000)
  hgs_out <- m_out(name = "Table S_FID_g - HU a", dep = "Escape distance", model = hgs, nsim = 5000)
  hgx_out <- m_out(name = "Table S_FID_g - HU b", dep = "Escape distance", model = hgx, nsim = 5000)
  ags_out <- m_out(name = "Table S_FID_g - AU a", dep = "Escape distance", model = ags, nsim = 5000)
  agx_out <- m_out(name = "Table S_FID_g - AU b", dep = "Escape distancey", model = agx, nsim = 5000)
  pgs_out <- m_out(name = "Table S_FID_g - PL a", dep = "Escape distance", model = pgs, nsim = 5000)
  pgx_out <- m_out(name = "Table S_FID_g - PL b", dep = "Escape distancey", model = pgx, nsim = 5000)

  out_FID_g <- rbind(mgs_out, mgx_out, fgs_out, fgx_out, pgs_out, pgx_out, cgs_out, cgx_out, hgs_out, hgx_out, ags_out, agx_out, fill = TRUE)
  out_FID_g[is.na(out_FID_g)] <- ""
  out_FID_g[, effect := gsub("scale\\(parks_percent_change_from_baseline\\)", "Google Mobility", effect)]
  out_FID_g[, effect := gsub("scale\\(Year\\)", "year", effect)]
  out_FID_g[, effect := gsub("scale\\(log\\(SD\\)", "starting distance (ln)", effect)]
  out_FID_g[, effect := gsub("scale\\(Temp\\)", "temperaturre", effect)]
  out_FID_g[, effect := gsub("scale\\(log\\(FlockSize\\)\\)", "flock size (ln)", effect)]
  out_FID_g[, effect := gsub("scale\\(log\\(BodyMass\\)\\)", "body mass (ln)", effect)]
  out_FID_g[, effect := gsub("scale\\(sin\\(rad\\)\\)", "sine of radians", effect)]
  out_FID_g[, effect := gsub("scale\\(cos\\(rad\\)\\)", "cosine of radians", effect)]
  out_FID_g[, effect := gsub("Species", "species", effect)]
  out_FID_g[, effect := gsub("sp_day_year", "species within day & year", effect)]
  out_FID_g[, effect := gsub("IDLocality", "site", effect)]
  out_FID_g[, effect := gsub("sp_loc", "species within site", effect)]
  fwrite(file = here::here("Outputs/Table_S_FID_g_c.csv"), out_FID_g)
load(here::here("Data/dat_est_Google_rev.Rdata"))
og[predictor %in% c("scale(parks_percent_change_from_baseline)"), predictor := "Google Mobility"]
ogo <- og[predictor %in% c("Google Mobility")]
ogo[, N:=as.numeric(sub('.*N = ', '', model))]
# add meta-analytical mean
  ogo_s = ogo[control_for_starting_distance == 'yes']
  met = summary(meta.summaries(d = ogo_s$estimate, se = ogo_s$sd, method = "fixed", weights = ogo_s$N))$summci
  ogo_met = data.table(predictor = "Period", estimate = met[2], lwr = met[1], upr = met[3], sd = NA, model = NA, control_for_starting_distance = "yes", Country = "Combined\n(metanalytical)", N = NA)

  ogo_sx = ogo[control_for_starting_distance == "no"]
  metx = summary(meta.summaries(d = ogo_sx$estimate, se = ogo_sx$sd, method = "fixed", weights = ogo_sx$N))$summci
  ogo_metx = data.table(predictor = "Period", estimate = metx[2], lwr = metx[1], upr = metx[3], sd = NA, model = NA, control_for_starting_distance = "no", Country = "Combined\n(metanalytical)", N = NA)
  
  ogo = rbind(ogo, ogo_met, ogo_metx)
    
ogo[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia", "Combined\n(metanalytical)", "All\n(mixed model)")))]

width_ <- .5 # spacing between error bars

#col_ <- c(brewer.pal(n = 12, name = "Paired"), "grey30", "grey80")
#Tol_bright <- c("#EE6677", "#228833", "#4477AA", "#CCBB44", "#66CCEE", "#AA3377", "#BBBBBB")
#Tol_muted <- c("#88CCEE", "#44AA99", "#117733", "#332288", "#DDCC77", "#999933", "#CC6677", "#882255", "#AA4499", "#DDDDDD")
#Tol_light <- c("#BBCC33", "#AAAA00", "#77AADD", "#EE8866", "#EEDD88", "#FFAABB", "#99DDFF", "#44BB99", "#DDDDDD")

# From Color Universal Design (CUD): https://jfly.uni-koeln.de/color/
#Okabe_Ito <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#000000")
#col_ = Okabe_Ito[7:1]
# JAMA and LocusZoom modified order
#col_ =  c("#374E55FF", "#374E55FF", "#DF8F44FF", "#79AF97FF", "#00A1D5FF", "#B24745FF",  "#80796BFF") #"#6A6599FF",
#col_ <- c("#357EBDFF", "#9632B8FF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#D43F3AFF", "#D43F3AFF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
col_ = c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
#show_col(col_)

#p = 
ggplot(ogo, aes(x = estimate, y = Country, col = Country, shape = control_for_starting_distance)) +
    geom_vline(xintercept = 0, color = "grey", linetype = "dotted") +
    geom_errorbarh(aes(xmin = lwr, xmax = upr), height = 0, position = ggstance::position_dodgev(width_)) +
    # geom_point(position = ggstance::position_dodgev(.6)) +
    geom_point(position = position_dodge(width = width_), bg = "white", size = 1.1) +
    # scale_color_viridis(discrete=TRUE, begin=0, end = 0.5)  +
    # scale_fill_viridis(discrete=TRUE, begin=0, end = 0.5) +
    # geom_text( aes(x = n_pos,label = N), vjust = 0, size = 1.75, position = ggstance::position_dodgev(width_))+ # 3 positions for 3 bars
    # annotate("text", x=log10(3), y=85, label= "Used", col = "grey30", size = 2.5)+

    scale_shape_manual(name = "Controlled for\nstarting distance", guide = guide_legend(reverse = TRUE), values = c(21, 19)) +
    #scale_color_jama(guide = "none")+ #, palette = 'light'
    scale_color_manual(guide = "none", values = col_) + #guide_legend(reverse = TRUE)
    scale_x_continuous(breaks = round(seq(-0.3, 0.2, by = 0.1), 1)) +
    ylab("") +
    xlab("Standardized effect size of Google Mobility (human presence)\n[on flight initiation distance]") +
    # coord_cartesian(xlim = c(-.15, .15)) +
    # scale_x_continuous(breaks = round(seq(-.15, .15, by = 0.05),2)) +
    theme_bw() +
    theme(
        legend.position = "right",
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 6),
        # legend.spacing.y = unit(0.1, 'cm'),
        legend.key.height = unit(0.5, "line"),
        legend.margin = margin(0, 0, 0, 0),
        # legend.position=c(0.5,1.6),
        plot.title = element_text(color = "grey", size = 7),
        plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 0.5, unit = "pt"),
        panel.grid = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = ax_lines, size = 0.25),
        axis.line.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
        # axis.text.x = element_text()
        axis.ticks.length = unit(1, "pt"),
        axis.text.x = element_text(, size = 6),
        axis.text.y = element_text(colour = "black", size = 7),
        axis.title = element_text(size = 7)
    )

ggsave(here::here("Outputs/Fig_Google_rev_width_CustomLocusZoom_v2.png"), width = 8, height = 6, unit = "cm", dpi = 600)

Fig. S_Google | The effect size with 95%CIs of Google Mobility on flight initiation distance (ln-transformed). The general model on all data was controlled for year, starting distance (ln-transformed; filled circles) or not (empty circles), flock size (ln-transformed), temperature (also a proxy for a day within the breeding season: r Pearson = 0.48; Fig. S1) and time of day. To account for circular properties of time, time was transformed into radians (2 × time × π/24) and fitted as sine and cosine of radians (Bulla et al. 2016). All continuous variables were standardised by subtracting the mean and dividing by the standard deviation. The multicollinearity was small as correlations between predictors were weak (r<XX, Fig. S1). To account for the non-independence of data points (Schielzeth & Forstmeier 2009; Barr et al. 2013), we fitted random intercepts of genus, species, species at a given day and year, site, and species within a site and in the full model also random slope of Google Mobility within genus and country. Note that in the full model, fitting Google Mobility as random slope at other random intercepts produces similar results (see Fig Sxx). We used this approach with a full dataset with all observations (n = xx), as well as with conservative datasets, one with at least five observations per species , the other with at least 10 observations per species. Albeit random structure of the full model accounts for the potential country specific effect (country accounted for xx% of variance, Table xx), depicted are also estimates from country specific models, with same predictors and random structure as the full model (excluding the random slopes and also random intercepts of country and in case of Poland site, as specific sites were not noted in Poland). Using meta.summaries function from rmeta R-package we used the country estimates, their standard deviation, and sample size per country to estimate the meta-analytical mean, which reflects the results based on the full mixed effect model. Importantly, the models test for within year changes in esccape distancces due to changes in human presence, i.e. day to day changes in Google Mobility. In other words, the models investigate whether birds flexibly adjust their fear response on a day to day basis (regardless of COVID-19 restrictions), a test different to the one intended in this study. The effect sizes are small, mostly close to zero, and hihgly uncertain.

Table FID_g | Escape distance in relations to Google Mobilitty

out_FID_g$R2_mar=out_FID_g$R2_con = NULL
out_FID_g %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
model response error_structure N type effect estimate_r lwr_r upr_r
Table S_FID_g - full a Escape distance Gaussian 3644 fixed (Intercept) 0.021 -0.095 0.142
fixed year 0.017 -0.024 0.059
fixed starting distance (ln)) 0.498 0.466 0.53
fixed flock size (ln) -0.002 -0.028 0.024
fixed body mass (ln) 0.027 -0.039 0.092
fixed sine of radians -0.009 -0.049 0.032
fixed cosine of radians 0.011 -0.022 0.043
fixed temperaturre 0.022 -0.015 0.06
fixed Google Mobility -0.096 -0.228 0.034
random % species within day & year (Intercept) 6% 5% 7%
random % species within site (Intercept) 5% 4% 5%
random % species (Intercept) 7% 7% 7%
random % site (Intercept) 9% 9% 9%
random % genus (Intercept) 0% -1% 4%
random % genus Google Mobility 0% -1% 4%
random % Country (Intercept) 1% -4% 4%
random % Country Google Mobility 1% -4% 4%
random % Residual 71% 59% 82%
Table S_FID_g - full b Escape distance Gaussian 3644 fixed (Intercept) 0.108 -0.143 0.361
fixed year 0.079 0.027 0.131
fixed flock size (ln) 0.012 -0.017 0.042
fixed body mass (ln) 0.132 0.049 0.217
fixed sine of radians 0.066 0.02 0.113
fixed cosine of radians 0.075 0.037 0.113
fixed temperaturre 0.053 0.005 0.1
fixed scale(StringencyIndex) 0.045 -0.028 0.118
fixed Google Mobility -0.117 -0.302 0.074
random % species within day & year (Intercept) 7% 5% 9%
random % species within site (Intercept) 6% 4% 8%
random % species (Intercept) 10% 8% 12%
random % site (Intercept) 13% 10% 15%
random % genus (Intercept) 0% -1% 5%
random % genus Google Mobility 0% -1% 5%
random % Country (Intercept) 1% -16% 10%
random % Country Google Mobility 1% -16% 10%
random % Residual 62% 41% 89%
Table S_FID_g - FI a Escape distance Gaussian 322 fixed (Intercept) 0.023 -0.194 0.232
fixed year 0.13 -0.011 0.272
fixed starting distance (ln)) 0.194 0.087 0.301
fixed flock size (ln) -0.119 -0.219 -0.015
fixed body mass (ln) 0.238 0.086 0.399
fixed sine of radians 0.132 -0.012 0.272
fixed cosine of radians 0.042 -0.081 0.163
fixed temperaturre -0.004 -0.139 0.131
fixed Google Mobility -0.1 -0.247 0.036
random % species within site (Intercept) 13% 12% 13%
random % species within day & year (Intercept) 3% 2% 3%
random % site (Intercept) 17% 13% 19%
random % species (Intercept) 0% 0% 0%
random % genus (Intercept) 2% 1% 4%
random % Residual 66% 62% 71%
Table S_FID_g - FI b Escape distance Gaussian 322 fixed (Intercept) -0.004 -0.226 0.225
fixed year 0.153 0.009 0.295
fixed flock size (ln) -0.113 -0.216 -0.009
fixed body mass (ln) 0.317 0.158 0.479
fixed sine of radians 0.119 -0.024 0.263
fixed cosine of radians 0.03 -0.091 0.156
fixed temperaturre -0.017 -0.153 0.116
fixed Google Mobility -0.088 -0.236 0.054
random % species within site (Intercept) 14% 14% 15%
random % species within day & year (Intercept) 3% 3% 3%
random % site (Intercept) 18% 15% 21%
random % species (Intercept) 0% 0% 1%
random % genus (Intercept) 2% 1% 4%
random % Residual 62% 57% 67%
Table S_FID_g - PL a Escape distance Gaussian 329 fixed (Intercept) 0.002 -0.127 0.135
fixed year -0.189 -0.327 -0.048
fixed starting distance (ln)) 0.532 0.446 0.62
fixed flock size (ln) -0.023 -0.107 0.061
fixed body mass (ln) -0.097 -0.203 0.013
fixed sine of radians 0.109 -0.14 0.357
fixed cosine of radians -0.163 -0.414 0.091
fixed temperaturre -0.175 -0.273 -0.077
fixed Google Mobility 0.004 -0.117 0.125
random % species within day & year (Intercept) 13% 13% 14%
random % species (Intercept) 7% 5% 10%
random % genus (Intercept) 0% 0% 0%
random % Residual 79% 76% 82%
Table S_FID_g - PL b Escape distancey Gaussian 329 fixed (Intercept) 0.001 -0.202 0.204
fixed year -0.378 -0.534 -0.222
fixed flock size (ln) 0.059 -0.037 0.158
fixed body mass (ln) -0.023 -0.178 0.126
fixed sine of radians -0.009 -0.31 0.284
fixed cosine of radians -0.036 -0.334 0.262
fixed temperaturre -0.159 -0.272 -0.045
fixed Google Mobility 0.054 -0.087 0.197
random % species within day & year (Intercept) 7% 7% 7%
random % species (Intercept) 21% 16% 26%
random % genus (Intercept) 0% 0% 0%
random % Residual 72% 67% 77%
Table S_FID_g - CZ a Escape distance Gaussian 1054 fixed (Intercept) 0.227 0.02 0.438
fixed starting distance (ln)) 0.345 0.287 0.402
fixed flock size (ln) 0.01 -0.042 0.059
fixed body mass (ln) 0.187 0.029 0.341
fixed sine of radians 0.062 -0.006 0.132
fixed cosine of radians 0.034 -0.022 0.089
fixed temperaturre 0.065 -0.006 0.135
fixed Google Mobility -0.089 -0.164 -0.016
random % species within site (Intercept) 5% 5% 5%
random % species within day & year (Intercept) 1% 0% 1%
random % species (Intercept) 14% 11% 17%
random % genus (Intercept) 14% 10% 18%
random % site (Intercept) 5% 4% 6%
random % Residual 61% 54% 69%
Table S_FID_g - CZ b Escape distance Gaussian 1054 fixed (Intercept) 0.229 0.025 0.436
fixed starting distance (ln)) 0.344 0.284 0.404
fixed flock size (ln) 0.009 -0.04 0.059
fixed body mass (ln) 0.182 0.031 0.343
fixed sine of radians 0.062 -0.009 0.131
fixed cosine of radians 0.035 -0.021 0.09
fixed temperaturre 0.065 -0.002 0.136
fixed Google Mobility -0.091 -0.163 -0.018
random % species within site (Intercept) 5% 5% 5%
random % species within day & year (Intercept) 1% 0% 1%
random % species (Intercept) 14% 12% 17%
random % genus (Intercept) 14% 11% 18%
random % site (Intercept) 5% 4% 6%
random % Residual 61% 54% 69%
Table S_FID_g - HU a Escape distance Gaussian 874 fixed (Intercept) 0.085 -0.183 0.35
fixed year 0.19 0.078 0.302
fixed starting distance (ln)) 0.506 0.443 0.567
fixed flock size (ln) 0.001 -0.053 0.053
fixed body mass (ln) 0 -0.142 0.139
fixed sine of radians -0.122 -0.191 -0.057
fixed cosine of radians 0.024 -0.041 0.088
fixed temperaturre -0.029 -0.109 0.05
fixed Google Mobility -0.028 -0.128 0.07
random % species within day & year (Intercept) 3% 3% 3%
random % species within site (Intercept) 0% 0% 0%
random % species (Intercept) 1% 1% 1%
random % genus (Intercept) 11% 7% 15%
random % site (Intercept) 9% 6% 13%
random % Residual 75% 67% 82%
Table S_FID_g - HU b Escape distance Gaussian 874 fixed (Intercept) 0.092 -0.353 0.52
fixed year 0.321 0.187 0.457
fixed flock size (ln) 0 -0.058 0.059
fixed body mass (ln) 0.284 0.066 0.511
fixed sine of radians -0.157 -0.233 -0.08
fixed cosine of radians 0.033 -0.039 0.104
fixed temperaturre -0.046 -0.138 0.052
fixed Google Mobility -0.107 -0.23 0.014
random % species within day & year (Intercept) 5% 5% 5%
random % species within site (Intercept) 0% 0% 0%
random % species (Intercept) 1% 1% 2%
random % genus (Intercept) 24% 18% 30%
random % site (Intercept) 16% 13% 18%
random % Residual 54% 45% 63%
Table S_FID_g - AU a Escape distance Gaussian 1065 fixed (Intercept) -0.056 -0.197 0.086
fixed year 0.028 -0.034 0.092
fixed starting distance (ln)) 0.499 0.448 0.552
fixed flock size (ln) 0.025 -0.024 0.074
fixed body mass (ln) -0.05 -0.145 0.044
fixed sine of radians -0.044 -0.117 0.028
fixed cosine of radians -0.043 -0.112 0.028
fixed temperaturre 0.019 -0.038 0.076
fixed Google Mobility -0.035 -0.092 0.021
random % species within day & year (Intercept) 6% 6% 6%
random % species within site (Intercept) 12% 12% 12%
random % species (Intercept) 12% 10% 13%
random % genus (Intercept) 2% 1% 2%
random % site (Intercept) 8% 6% 10%
random % Residual 61% 57% 65%
Table S_FID_g - AU b Escape distancey Gaussian 1065 fixed (Intercept) -0.108 -0.302 0.08
fixed year 0.082 0.01 0.157
fixed flock size (ln) 0.055 0 0.11
fixed body mass (ln) 0.068 -0.059 0.196
fixed sine of radians -0.031 -0.116 0.054
fixed cosine of radians -0.017 -0.097 0.063
fixed temperaturre 0.007 -0.056 0.073
fixed Google Mobility -0.018 -0.082 0.045
random % species within day & year (Intercept) 7% 7% 7%
random % species within site (Intercept) 11% 10% 11%
random % species (Intercept) 14% 12% 16%
random % genus (Intercept) 6% 5% 7%
random % site (Intercept) 10% 8% 12%
random % Residual 52% 48% 58%

Continuous variables were z-transformed.

Genus and spcies ‘raw’ data visualisations

dxx <- d[paste(IDLocality, Species) %in% paste(pp$IDLocality, pp$Species)]
# length(dxx[, unique(paste(sp_loc, Covid))])
m <- lm(log(FID) ~ log(SD), dxx)
dxx[, resid_FID := resid(m)]
a <- dxx[, .(mean(resid_FID), sd(resid_FID), mean(FID), .N), by = .(Country, IDLocality, genus, Species, sp_loc, Covid)]
setnames(a, old = c("V1", "V2", "V3"), new = c("resid_FID_avg", "SD", "FID_avg"))
a[is.na(SD), SD := 0]

aw <- reshape(a, idvar = c("Country", "IDLocality", "genus", "Species", "sp_loc"), timevar = "Covid", direction = "wide")
aw[, Species := gsub("[_]", " ", Species)]
aw <- merge(aw, t, all.x = TRUE)
#table(aw$Family)

x <- aw[, .N, by = Species]
#x[order(Species)]
aw[, genus2 := genus]
aw[Species %in% x[N %in% c(1, 2), Species], genus2 := "other"]

#aw[genus2 == "Phoenicurus", unique(Species)]

ph[genus2 == "Motacilla" | uid %in% c("67a9ecfd-58ba-44a4-9986-243b6e610419"), uid := "cf522e02-35cc-44f5-841c-0e642987c2e4"]
ph[genus2 == "Sylvia", uid := "67a9ecfd-58ba-44a4-9986-243b6e610419"]

ph[, size := 0.2]
ph[genus2 %in% c("Anas", "Columba", "Dendrocopos", "Sturnus"), size := c(0.25, 0.25, 0.15, 0.1)]
ph[, FID_avg.0 := 1.5]
ph[, FID_avg.1 := 20]
ph[genus2 %in% c("Anas", "Columba"), FID_avg.0 := c(1.7, 1.7)]

ph[, resid_FID_avg.0 := -1.7]
ph[, resid_FID_avg.1 := 0.7]
ph[genus2 %in% c("Anas", "Columba"), resid_FID_avg.0 := c(-1.6, -1.6)]

ph[, genus2 := factor(genus2, levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other"))]

aw[, genus2 := factor(genus2, levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other"))]

aw[, genus2 := factor(genus2, levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other"))]

# Fig?? and left panel of S5- plot from files
anas <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Anas.png"))))
columba <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Columba.png"))))
Dendrocopos <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Dendrocopos.png"))))
Larus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Larus_flip.png"))))
Picus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Picus.png"))))
Motacilla <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Motacilla.png"))))
Erithacus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Erithacus.png"))))
Phoenicurus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Phoenicurus.png"))))
Turdus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Turdus.png"))))
Sylvia <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Sylvia.png"))))
Parus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Parus_flip.png"))))
Sitta <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Sitta.png"))))
Pica <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Pica.png"))))
Garrulus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Garrulus.png"))))
Corvus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Corvus.png"))))
Sturnus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Sturnus.png"))))
Passer <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Passer_flip.png"))))
Fringilla <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Fringilla.png"))))
other <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/other_flip.png"))))

ann_text <- data.frame(
  FID_avg.0 = 8, FID_avg.1 = 10, lab = "Text",
  genus2 = factor("Anas", levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other"))
)
ann_text2 <- data.frame(
  FID_avg.0 = 6, FID_avg.1 = 3, lab = "Text",
  genus2 = factor("Larus", levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other"))
)

aw2 <- data.frame(FID_avg.0 = c(11.25, 11.25), FID_avg.1 = c(3.5, 5.8), genus2 = factor("Larus", levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other")))

col3_ <- c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1]
col3__ <- col3_[3:7]

aw[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia")))]

g_gen <-
  ggplot(aw, aes(x = FID_avg.0, y = FID_avg.1)) +
  # geom_errorbar(aes(ymin = FID_avg.1-SD.1, ymax = FID_avg.1+SD.1, col = Country), width = 0) +
  # geom_errorbar(aes(xmin = FID_avg.0-SD.0, xmax = FID_avg.0+SD.0, col = Country), width = 0) +
  # geom_point(pch = 21, alpha = 0.7, aes(col = Country)) +
  annotation_custom2(anas, data = ph[genus2 == "Anas"], xmin = 0.05, xmax = 0.5, ymax = 2.6) +
  annotation_custom2(Larus, data = ph[genus2 == "Larus"], xmin = 0.05, xmax = 0.5, ymax = 2.6) +
  annotation_custom2(columba, data = ph[genus2 == "Columba"], xmin = 0.05, xmax = 0.4, ymax = 2.7) +
  annotation_custom2(Dendrocopos, data = ph[genus2 == "Dendrocopos"], xmin = 0.05, xmax = 0.25, ymax = 2.6) +
  annotation_custom2(Picus, data = ph[genus2 == "Picus"], xmin = 0.05, xmax = 0.4, ymax = 2.7) +
  annotation_custom2(Motacilla, data = ph[genus2 == "Motacilla"], xmin = 0.05, xmax = 0.5, ymax = 2.7) +
  annotation_custom2(Erithacus, data = ph[genus2 == "Erithacus"], xmin = 0.05, xmax = 0.35, ymax = 2.7) +
  annotation_custom2(Phoenicurus, data = ph[genus2 == "Phoenicurus"], xmin = 0.05, xmax = 0.35, ymax = 2.7) +
  annotation_custom2(Turdus, data = ph[genus2 == "Turdus"], xmin = 0.05, xmax = 0.5, ymax = 2.7) +
  annotation_custom2(Sylvia, data = ph[genus2 == "Sylvia"], xmin = 0.05, xmax = 0.5, ymax = 2.7) +
  annotation_custom2(Parus, data = ph[genus2 == "Parus"], xmin = 0.05, xmax = 0.42, ymax = 2.7) +
  annotation_custom2(Sitta, data = ph[genus2 == "Sitta"], xmin = 0.05, xmax = 0.5, ymax = 2.7) +
  annotation_custom2(Pica, data = ph[genus2 == "Pica"], xmin = 0.05, xmax = 0.5, ymax = 2.5) +
  annotation_custom2(Garrulus, data = ph[genus2 == "Garrulus"], xmin = 0.05, xmax = 0.6, ymax = 2.7) +
  annotation_custom2(Corvus, data = ph[genus2 == "Corvus"], xmin = 0.05, xmax = 0.4, ymax = 2.55) +
  annotation_custom2(Sturnus, data = ph[genus2 == "Sturnus"], xmin = 0.05, xmax = 0.24, ymax = 2.65) +
  annotation_custom2(Passer, data = ph[genus2 == "Passer"], xmin = 0.05, xmax = 0.36, ymax = 2.65) +
  annotation_custom2(Fringilla, data = ph[genus2 == "Fringilla"], xmin = 0.05, xmax = 0.5, ymax = 2.7) +
  annotation_custom2(other, data = ph[genus2 == "other"], xmin = 0.05, xmax = 0.45, ymax = 2.4) +
  geom_point(pch = 21, alpha = 0.7, aes(fill = Country), col = 'grey20', alpha =0.8)+#col = "white") +
  # ggtitle ("Sim based")+
  geom_abline(intercept = 0, slope = 1, lty = 3, col = "grey80") +
  geom_line(data = aw2, col = "grey80", lwd = 0.25) +
  geom_text(data = ann_text, label = "No difference", col = "grey80", angle = 45, size = 2) +
  geom_text(data = ann_text2, label = "Species mean / site", col = "grey60", size = 2, ) +
  facet_wrap(~genus2) +
  # geom_phylopic(data = o, aes(image = uid),  color = "grey80", size = o$size) + # ,
  #scale_fill_viridis(discrete = TRUE, guide = guide_legend(reverse = FALSE)) +
  scale_fill_manual(values = col3__,guide = guide_legend(reverse = TRUE)) +
  scale_x_continuous("Before COVID-19 shutdown - flight initiation distance [m]", expand = c(0, 0), trans = "log10") +
  scale_y_continuous("During COVID-19 shutdown - flight initiation distance [m]", expand = c(0, 0), trans = "log10") +
  # labs(title = "Species means per sampling location")+
  theme_MB +
  theme(
    plot.title = element_text(size = 7),
    strip.background = element_blank(),
    # panel.spacing = unit(1, "mm"),
    legend.position = c(1, 0.025),
    legend.justification = c(1, -0.05)
  )
gg_gen <- ggplotGrob(g_gen) # gg$layout$name
ggx_gen <- gtable_filter_remove(gg_gen,
  name = paste0("axis-b-", c(2, 4), "-4"),
  trim = FALSE
)
# grid.draw(ggx)
# ggsave('Outputs/Fig_2_width-114mm.png',ggx, width=4.5,height=4.5,dpi=600) # 11.43cm # with label on top
if(save_plot==TRUE){
ggsave(here::here("Outputs/Fig_2_width-122mm_col_grey_rev.png"), ggx_gen, width = 4.8, height = 4.5, dpi = 600) # 12.2cm # with label inside
}
   # Fig S5 right panel
   g_2 <-
     ggplot(aw, aes(x = resid_FID_avg.0, y = resid_FID_avg.1)) +
     geom_point(pch = 21, alpha = 0.7, aes(fill = Country), col = 'grey20', alpha =0.8)+#col = "white") +
     geom_abline(intercept = 0, slope = 1, lty = 3, col = "grey80") +
     facet_wrap(~genus2) +
     # geom_phylopic(data = o, aes(image = uid),  color = "grey80", size = o$size) +
     #scale_fill_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
     scale_fill_manual(values = col3__,guide = guide_legend(reverse = TRUE)) +
     scale_x_continuous("Before COVID-19 shutdown - residual escape distance", expand = c(0, 0)) +
     scale_y_continuous("During COVID-19 shutdown - residual escape distance", expand = c(0, 0)) +
     # labs(title = "Species means per sampling location")+
     theme_MB +
     theme(
       plot.title = element_text(size = 7),
       strip.background = element_blank(),
       legend.position = "none",
       # legend.position = c(1, 0),
       legend.justification = c(1, 0)
     )

   g_g2 <- ggplotGrob(g_2) # gg$layout$name
   g_gx2 <- gtable_filter_remove(g_g2,
     name = paste0("axis-b-", c(2, 4), "-4"),
     trim = FALSE
   )
   # Fig S5 combine
   grid.draw(cbind(ggx_gen, g_gx2, size = "last"))

   if(save_plot==TRUE){
   ggsave(here::here("Outputs/Fig_S5_rev_v2.png"), cbind(ggx_gen, g_gx2, size = "last"), width = 4.8 * 2, height = 4.5, dpi = 600)
   }

exploration for Fig S5 legend

correlation between mean fid and residual fid

   #ggplot(a, aes(x = resid_FID_avg)) + geom_histogram()
   g_a = 
   ggplot(a, aes(x = FID_avg, y = resid_FID_avg)) +
     labs(subtitle='Species means per period')+
     stat_smooth() +
     geom_point() +
     stat_cor(aes(label = ..r.label..), label.x = 0.3, size = 2) +
     scale_x_continuous(trans = "log")
   g_r = 
   ggplot(dxx, aes(x = FID, y = resid_FID)) +
     labs(subtitle='Raw observations')+
     stat_smooth() +
     geom_point() +
     stat_cor(aes(label = ..r.label..), label.x = 0.3, size = 2) +
     scale_x_continuous(trans = "log")
   grid.draw(cbind(ggplotGrob(g_a), ggplotGrob(g_r)))

   #cor(log(dxx$FID), dxx$resid_FID)
   #cor(log(a$FID_avg), a$resid_FID_avg)

difference in number of observations per species and site before aand during shutdowns

   ggplot(aw, aes(x = N.0 - N.1)) +
     geom_histogram()

   #nrow(aw[abs(N.0 - N.1) > 2])
   #nrow(aw[!abs(N.0 - N.1) > 2])

Fig S3 - species means

      aw[, sp2 := gsub(" ", "\n", Species)]
      ann_text <- data.frame(FID_avg.0 = 8, FID_avg.1 = 10,lab = "Text",
                       Species = factor('Aegithalos caudatus',levels = levels(as.factor(aw$Species))))
      ann_text$sp2 = gsub(" ", "\n", ann_text$Species)
      g3 = 
        ggplot(aw, aes(x = FID_avg.0, y = FID_avg.1)) + 
          #geom_errorbar(aes(ymin = FID_avg.1-SD.1, ymax = FID_avg.1+SD.1, col = Country), width = 0) +
          #geom_errorbar(aes(xmin = FID_avg.0-SD.0, xmax = FID_avg.0+SD.0, col = Country), width = 0) +
          #geom_point(pch = 21, alpha = 0.7, aes(col = Country)) + 
          geom_point(pch = 21, alpha = 0.7, aes(fill = Country), col = 'grey20', alpha =0.8)+#col = "white") +
            #ggtitle ("Sim based")+
          geom_abline(intercept = 0, slope = 1, lty =3, col = "grey80")+
          geom_text(data = ann_text,label = "No difference", col = "grey80",angle = 45, size = 2) + 
          facet_wrap(~sp2) +
          #geom_phylopic(data = o, aes(image = uid),  color = "grey80", size = o$size) + # ,
          #scale_fill_viridis(discrete=TRUE,guide = guide_legend(reverse = FALSE))  +
          scale_fill_manual(values = col3__,guide = guide_legend(reverse = TRUE)) +
          scale_x_continuous("Before COVID-19 shutdown - flight initiation distance [m]", expand = c(0, 0), trans = 'log10') +
          scale_y_continuous("During COVID-19 shutdown - flight initiation distance [m]", expand = c(0, 0), trans = 'log10') +
          labs(title = "Species means per sampling location")+
          theme_MB  +
          theme(
                  plot.title = element_text(size=7),
                  strip.background = element_blank(),
                  #panel.spacing = unit(1, "mm"),
                  legend.position = c(0.96, 0.0),
                  legend.justification = c(1, 0)
                  )  
        gg3 <- ggplotGrob(g3) #gg2$layout$name
        ggx3 <- gtable_filter_remove(gg3, name = c(paste0("axis-b-", c(2, 4), "-7"), "axis-b-6-6"),
                                         trim = FALSE)
        grid.draw(ggx3)

        if(save_plot==TRUE){
        ggsave(here::here('Outputs/Fig_S3_species_rev_v2.png'),ggx3, width=13.5,height=17.5,unit = 'cm', dpi=600) # 11.43cm
        }

Figure S2 - Alternative models give similar results

  # prepare estimates Period
    # 01a all data, main text
      m1a <- lmer(scale(log(FID)) ~
        scale(log(SD)) +
        scale(log(FlockSize)) +
        scale(log(BodyMass)) +
        scale(sin(rad)) + scale(cos(rad)) +
        # scale(Day)+
        scale(Temp) +
        scale(Covid) +
        (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | Country) + (scale(Covid) | IDLocality) + (1 | sp_loc),
      data = d, REML = FALSE,
       control = lmerControl(
         optimizer = "optimx", optCtrl = list(method = "nlminb")
       )
      )
     est_m1a = est_out(m1a, "01a) (1|Year) + (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (1|sp_loc)")
    # 01b all data, all random slopes - singularity 
      m1b=lmer(scale(log(FID))~
                  scale(log(SD))+
                  scale(log(FlockSize))+
                  scale(log(BodyMass))+
                  scale(sin(rad)) + scale(cos(rad)) + 
                  #scale(Day)+
                  scale(Temp)+
                  scale(Covid)+
                  (1|Year) + (1 | weekday) + (scale(Covid)|genus)+(scale(Covid)|Species)+(1|sp_day_year) + (scale(Covid)|Country) +(scale(Covid)|IDLocality) + (scale(Covid)|sp_loc),
                  data = d, REML = FALSE
                )
                  # # (Covid|IDLocality) +
      est_m1b = est_out(m1b, "01b) (1|Year) + (1|weekday) + (scale(Covid)|genus) + (scale(Covid)|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc)")
    # 01c all data, all random slopes, but some without cor to avoid singularity
      m1c=lmer(scale(log(FID))~
                  scale(log(SD))+
                  scale(log(FlockSize))+
                  scale(log(BodyMass))+
                  scale(sin(rad)) + scale(cos(rad)) + 
                  #scale(Day)+
                  scale(Temp)+
                  scale(Covid)+
                  #(1|Year) +(0+Covid|genus)+(0+Covid|Species)+(1|sp_day_year) + (Covid|Country) + (0+Covid|IDLocality) +(Covid|sp_loc)
                  (1|Year)+ (1|weekday) + (0+scale(Covid)|genus)+(0+scale(Covid)|Species)+(1|sp_day_year) + (scale(Covid)|Country) +(scale(Covid)|IDLocality) + (scale(Covid)|sp_loc),
                  data = d, REML = FALSE) # (0+scale(Covid)|genus) is zero 
      est_m1c = est_out(m1c, "01c) (1|Year) + (1|weekday) + (0+scale(Covid)|genus) + (0+scale(Covid)|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc)")
    # 01d all data, random slopes that allow for non-singular fit 
      m1d=lmer(scale(log(FID))~
                  scale(log(SD))+
                  scale(log(FlockSize))+
                  scale(log(BodyMass))+
                  scale(sin(rad)) + scale(cos(rad)) + 
                  #scale(Day)+
                  scale(Temp)+
                  scale(Covid)+
                  (1|Year) + (1|weekday) +(1|genus)+(1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc),
                  data = d, REML = FALSE)# (Covid|IDLocality) +
      #d[,res := resid(mf)]
      est_m1d = est_out(m1d, "01d) (1|Year) + (1|weekday)+ (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc)")
    # 01e all data, random slopes that allow for non-singular fit with simple structure
      m1e=lmer(scale(log(FID))~
                  scale(log(SD))+
                  scale(log(FlockSize))+
                  scale(log(BodyMass))+
                  scale(sin(rad)) + scale(cos(rad)) + 
                  #scale(Day)+
                  scale(Temp)+
                  scale(Covid)+
                  (1|Year) + (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc),
                  data = d, REML = FALSE)# (Covid|IDLocality) +
      est_m1e = est_out(m1e, "01e) (1|Year) + (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc)")

    # 02a) before & during > 4/species - main text 
      dx = dd[N_during>4 & N_before >4]
      #dxx = d[Species %in% dx$Species]
      #dx2 = dd[N_during>9 & N_before >9]
      m2a <- lmer(scale(log(FID)) ~
        scale(log(SD)) +
        scale(log(FlockSize)) +
        scale(log(BodyMass)) +
        scale(sin(rad)) + scale(cos(rad)) +
        # scale(Day)+
        scale(Temp) +
        scale(Covid) +
        (1 | Year) + (1| weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | Country) + (scale(Covid) | IDLocality) + (1 | sp_loc),
      data = d[Species %in% dx$Species], REML = FALSE,
       control = lmerControl(
         optimizer = "optimx", optCtrl = list(method = "nlminb")
       )
      )
      est_m2a = est_out(m2a, "02a) (1|Year) + (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (1|sp_loc); >4/species/period")
    
    # 02b) before & during > 4/species - singularity
      dx = dd[N_during > 4 & N_before > 4]
      m2b=lmer(scale(log(FID))~
                  scale(log(SD))+
                  scale(log(FlockSize))+
                  scale(log(BodyMass))+
                  scale(sin(rad)) + scale(cos(rad)) + 
                  #scale(Day)+
                  scale(Temp)+
                  scale(Covid)+
                  (1|Year) + (1| weekday)+ (scale(Covid)|genus)+(scale(Covid)|Species)+(1|sp_day_year) + (scale(Covid)|Country) +(scale(Covid)|IDLocality) + (scale(Covid)|sp_loc),
                  #(1|Year) + (1|genus)+(1|Species)+(1|sp_day_year) + (Covid|Country) + (Covid|sp_loc),
                  data = d[Species %in% dx$Species],
                  REML = FALSE) # (Covid|IDLocality) +
      est_m2b = est_out(m2b, "02b) (1|Year) + (1|weekday)+ (scale(Covid)|genus)+(scale(Covid)|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc); >4/species/period")
    # 02c) before & during > 4/species, random slopes that allow for non-singular fit and simple structure
      dx = dd[N_during>4 & N_before >4]
      #dxx = d[Species %in% dx$Species]
      #dx2 = dd[N_during>9 & N_before >9]
      m2c=lmer(scale(log(FID))~
                  scale(log(SD))+
                  scale(log(FlockSize))+
                  scale(log(BodyMass))+
                  scale(sin(rad)) + scale(cos(rad)) + 
                  #scale(Day)+
                  scale(Temp)+
                  scale(Covid)+
                  (1|Year)+ (1| weekday) +(1|genus)+(1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc),
                  #(1|Year) + (1|genus)+(1|Species)+(1|sp_day_year) + (Covid|Country) + (Covid|sp_loc),
                  data = d[Species %in% dx$Species],
                  REML = FALSE) # (Covid|IDLocality) +
      est_m2c = est_out(m2c, "02c) (1|Year) + (1|weekday)+ (1|genus) + (1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc); >4/species/period")
    # 03a) before & during > 9/species - main text
      dx = dd[N_during > 9 & N_before > 9]
      m3a <- lmer(scale(log(FID)) ~
       scale(log(SD)) +
       scale(log(FlockSize)) +
       scale(log(BodyMass)) +
       scale(sin(rad)) + scale(cos(rad)) +
       # scale(Day)+
       scale(Temp) +
       scale(Covid) +
       (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | Country) + (scale(Covid) | IDLocality) + (1 | sp_loc),
      data = d[Species %in% dx$Species], REML = FALSE,
      )
     est_m3a = est_out(m3a, "03a) (1|Year) + (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (1|sp_loc); >9/species/period")
    # 03b) before & during > 9/species - singularity 
      dx = dd[N_during>9 & N_before >9]
      #dxx = d[Species %in% dx$Species]
      #dx2 = dd[N_during>9 & N_before >9]
      m3b=lmer(scale(log(FID))~
                  scale(log(SD))+
                  scale(log(FlockSize))+
                  scale(log(BodyMass))+
                  scale(sin(rad)) + scale(cos(rad)) + 
                  #scale(Day)+
                  scale(Temp)+
                  scale(Covid)+
                  (1|Year) + (1| weekday) + (scale(Covid)|genus)+(scale(Covid)|Species)+(1|sp_day_year) + (scale(Covid)|Country) +(scale(Covid)|IDLocality) + (scale(Covid)|sp_loc),
                  #(1|Year) + (1|genus)+(1|Species)+(1|sp_day_year) + (Covid|Country) + (Covid|sp_loc),
                  data = d[Species %in% dx$Species],
                  REML = FALSE) # (Covid|IDLocality) +
      est_m3b = est_out(m3b, "03b) (1|Year) + (1|weekday) + (scale(Covid)|genus)+(scale(Covid)|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc); >9/species/period")
    # 03c) before & during > 9/species, random slopes that allow for non-singular fit and simple structure
      dx = dd[N_during>9 & N_before >9]
      #dxx = d[Species %in% dx$Species]
      #dx2 = dd[N_during>9 & N_before >9]
      m3c=lmer(scale(log(FID))~
                  scale(log(SD))+
                  scale(log(FlockSize))+
                  scale(log(BodyMass))+
                  scale(sin(rad)) + scale(cos(rad)) + 
                  #scale(Day)+
                  scale(Temp)+
                  scale(Covid)+
                  (1|Year) + (1| weekday) +(1|genus)+(1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc),
                  #(1|Year) + (1|genus)+(1|Species)+(1|sp_day_year) + (Covid|Country) + (Covid|sp_loc),
                  data = d[Species %in% dx$Species],
                  REML = FALSE) # (Covid|IDLocality) +
      est_m3c = est_out(m3c, '03c) (1|Year) + (1|genus) + (1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc); >9/species/period')  
  
  # prepare estimates Stringency 
    m01a <- lmer(scale(log(FID)) ~
      scale(Year) +
      scale(log(SD)) +
      scale(log(FlockSize)) +
      scale(log(BodyMass)) +
      scale(sin(rad)) + scale(cos(rad)) +
      # scale(Day)+
      scale(Temp) +
      scale(StringencyIndex) +
      (1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1 | IDLocality) + (1 | sp_loc),
      data = s, REML = FALSE)
    est_m01a = est_out(m01a, "01a) (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc)")

    m01b=lmer(scale(log(FID))~ 
        scale(Year)+ 
        scale(log(SD))+ 
        scale(log(FlockSize))+ 
        scale(log(BodyMass))+ 
        scale(sin(rad)) + scale(cos(rad)) +  
        #scale(Day)+ 
        scale(Temp)+ 
        scale(StringencyIndex)+ 
        (1|weekday) + (scale(StringencyIndex)|genus)+(1|Species)+(1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc),
        data = s, REML = FALSE
        )  
        # (1|Year) explains nothing - could stay 
     est_m01b = est_out(m01b, "01b) (1|weekday) + (scale(StringencyIndex)|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc)")
     
     m01c=lmer(scale(log(FID))~
          scale(Year)+       
          scale(log(SD))+
          scale(log(FlockSize))+
          scale(log(BodyMass))+
          scale(sin(rad)) + scale(cos(rad)) + 
          #scale(Day)+
          scale(Temp)+
          scale(StringencyIndex)+
          (scale(StringencyIndex)|Country) + (1|IDLocality),  
          data = s, REML = FALSE) 
          # (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc),  explain nothing - could stay
     est_m01c = est_out(m01c, "01c) (scale(StringencyIndex)|Country) + (1|IDLocality)")

     m02a <- lmer(scale(log(FID)) ~
       scale(Year) +
       scale(log(SD)) +
       scale(log(FlockSize)) +
       scale(log(BodyMass)) +
       scale(sin(rad)) + scale(cos(rad)) +
       # scale(Day)+
       scale(Temp) +
       scale(StringencyIndex) +
       (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1 | IDLocality) + (1 | sp_loc),
     data = s[Nsp>4], REML = FALSE
     )

     est_m02a = est_out(m02a, "02a) (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc); >4/species")
    m02b=lmer(scale(log(FID))~ 
        scale(Year)+ 
        scale(log(SD))+ 
        scale(log(FlockSize))+ 
        scale(log(BodyMass))+ 
        scale(sin(rad)) + scale(cos(rad)) +  
        #scale(Day)+ 
        scale(Temp)+ 
        scale(StringencyIndex)+ 
        (1|weekday) + (scale(StringencyIndex)|genus)+(1|Species)+(1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc),
        data = s[Nsp>4], REML = FALSE,  
        control = lmerControl( 
            optimizer ='optimx', optCtrl=list(method='nlminb')) 
        )  
        # (1|Year) explains nothing - could stay 
     est_m02b = est_out(m02b, "02b) (1|weekday) + (scale(StringencyIndex)|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc); >4/species")
     
     m02c=lmer(scale(log(FID))~
          scale(Year)+       
          scale(log(SD))+
          scale(log(FlockSize))+
          scale(log(BodyMass))+
          scale(sin(rad)) + scale(cos(rad)) + 
          #scale(Day)+
          scale(Temp)+
          scale(StringencyIndex)+
          (scale(StringencyIndex)|Country) + (1|IDLocality),  
          data = s[Nsp>4], REML = FALSE) 
          # (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc),  explain nothing - could stay
      est_m02c = est_out(m02c, "02c) (scale(StringencyIndex)|Country) + (1|IDLocality); >4/species")

      m03a <- lmer(scale(log(FID)) ~
           scale(Year) +
           scale(log(SD)) +
           scale(log(FlockSize)) +
           scale(log(BodyMass)) +
           scale(sin(rad)) + scale(cos(rad)) +
           # scale(Day)+
           scale(Temp) +
           scale(StringencyIndex) +
           (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1 | IDLocality) + (1 | sp_loc),
         data = s[Nsp > 9], REML = FALSE
         )

         est_m03a = est_out(m03a, "03a) (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc); >9/species")
         m03b = lmer(scale(log(FID)) ~
           scale(Year) +
           scale(log(SD)) +
           scale(log(FlockSize)) +
           scale(log(BodyMass)) +
           scale(sin(rad)) + scale(cos(rad)) +
           # scale(Day)+
           scale(Temp) +
           scale(StringencyIndex) +
           (1 | weekday) + (scale(StringencyIndex) | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1 | IDLocality) + (1 | sp_loc),
         data = s[Nsp > 9], REML = FALSE
         )
         est_m03b = est_out(m03b, "03b) (1|weekday) + (scale(StringencyIndex)|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc); >9/species")

        m03c = lmer(scale(log(FID)) ~
           scale(Year) +
           scale(log(SD)) +
           scale(log(FlockSize)) +
           scale(log(BodyMass)) +
           scale(sin(rad)) + scale(cos(rad)) +
           # scale(Day)+
           scale(Temp) +
           scale(StringencyIndex) +
           (scale(StringencyIndex) | Country) + (1 | IDLocality),
         data = s[Nsp > 9], REML = FALSE
         )
         # (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc),  explain nothing - could stay
         est_m03c = est_out(m03c, "03c) (scale(StringencyIndex)|Country) + 1|IDLocality); >9/species")

  # prepare estimates google mobility 
    mg01a=lmer(scale(log(FID))~
      scale(Year)+       
      scale(log(SD))+
      scale(log(FlockSize))+
      scale(log(BodyMass))+
      scale(sin(rad)) + scale(cos(rad)) + 
      #scale(Day)+
      scale(Temp)+
      scale(parks_percent_change_from_baseline)+
      (1|weekday) + (1|genus) +(1|Species) + (1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality)+(1|sp_loc),  
      data = ss, REML = FALSE
      ) 
      # (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc),  explain nothing - could stay
    est_mg01a = est_out(mg01a, "01a) (1|weekday) + (1|genus) +(1|Species) + (1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality)+(1|sp_loc)")
    
    mg01b <- lmer(scale(log(FID)) ~
     scale(Year) +
     scale(log(SD)) +
     scale(log(FlockSize)) +
     scale(log(BodyMass)) +
     scale(sin(rad)) + scale(cos(rad)) +
     # scale(Day)+
     scale(Temp) +
     scale(parks_percent_change_from_baseline) +
     (1|weekday) + (scale(parks_percent_change_from_baseline)| genus) + (1 | Species) + (1 | sp_day_year) + 
     (scale(parks_percent_change_from_baseline)| Country) + (1 | IDLocality) + (1 | sp_loc),
     data = ss, REML = FALSE
     ) 

     est_mg01b = est_out(mg01b, "01b) (1|weeekday) + (scale(parks_percent_change_from_baseline)|genus)+(1|Species)+(1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality) +(1|sp_loc)")
     
    mg01c <- lmer(scale(log(FID)) ~
      scale(Year) +
      scale(log(SD)) +
      scale(log(FlockSize)) +
      scale(log(BodyMass)) +
      scale(sin(rad)) + scale(cos(rad)) +
      # scale(Day)+
      scale(Temp) +
      scale(parks_percent_change_from_baseline) +
      (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality),
    data = ss, REML = FALSE
    )
    est_mg01c = est_out(mg01c, "01c) (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality)")

    mg02a = lmer(scale(log(FID)) ~
      scale(Year) +
      scale(log(SD)) +
      scale(log(FlockSize)) +
      scale(log(BodyMass)) +
      scale(sin(rad)) + scale(cos(rad)) +
      # scale(Day)+
      scale(Temp) +
      scale(parks_percent_change_from_baseline) +
      (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality) + (1 | sp_loc),
    data = ss[Nsp>4], REML = FALSE
    )
    # (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc),  explain nothing - could stay
    est_mg02a = est_out(mg02a, "02a) (1|weekday) + (1|genus) +(1|Species) + (1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality)+(1|sp_loc); >4/specie")

    mg02b <- lmer(scale(log(FID)) ~
      scale(Year) +
      scale(log(SD)) +
      scale(log(FlockSize)) +
      scale(log(BodyMass)) +
      scale(sin(rad)) + scale(cos(rad)) +
      # scale(Day)+
      scale(Temp) +
      scale(parks_percent_change_from_baseline) +
      (1 | weekday) + (scale(parks_percent_change_from_baseline) | genus) + (1 | Species) + (1 | sp_day_year) +
      (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality) + (1 | sp_loc),
    data = ss[Nsp>4], REML = FALSE
    )

    est_mg02b = est_out(mg02b, "02b) (1|weeekday) + (scale(parks_percent_change_from_baseline)|genus)+(1|Species)+(1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality) +(1|sp_loc); >4/specie")

    mg02c <- lmer(scale(log(FID)) ~
      scale(Year) +
      scale(log(SD)) +
      scale(log(FlockSize)) +
      scale(log(BodyMass)) +
      scale(sin(rad)) + scale(cos(rad)) +
      # scale(Day)+
      scale(Temp) +
      scale(parks_percent_change_from_baseline) +
      (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality),
    data = ss[Nsp>4], REML = FALSE
    )
    est_mg02c = est_out(mg02c, "02c) (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality); >4/specie")

    mg03a = lmer(scale(log(FID)) ~
      scale(Year) +
      scale(log(SD)) +
      scale(log(FlockSize)) +
      scale(log(BodyMass)) +
      scale(sin(rad)) + scale(cos(rad)) +
      # scale(Day)+
      scale(Temp) +
      scale(parks_percent_change_from_baseline) +
      (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality) + (1 | sp_loc),
    data = ss[Nsp > 9], REML = FALSE
    )
    # (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc),  explain nothing - could stay
    est_mg03a = est_out(mg03a, "03a) (1|weekday) + (1|genus) +(1|Species) + (1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality)+(1|sp_loc); >9/specie")

    mg03b <- lmer(scale(log(FID)) ~
      scale(Year) +
      scale(log(SD)) +
      scale(log(FlockSize)) +
      scale(log(BodyMass)) +
      scale(sin(rad)) + scale(cos(rad)) +
      # scale(Day)+
      scale(Temp) +
      scale(parks_percent_change_from_baseline) +
      (1 | weekday) + (scale(parks_percent_change_from_baseline) | genus) + (1 | Species) + (1 | sp_day_year) +
      (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality) + (1 | sp_loc),
    data = ss[Nsp > 9], REML = FALSE
    )

    est_mg03b = est_out(mg03b, "03b) (1|weeekday) + (scale(parks_percent_change_from_baseline)|genus)+(1|Species)+(1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality) +(1|sp_loc); >9/specie")

    mg03c <- lmer(scale(log(FID)) ~
      scale(Year) +
      scale(log(SD)) +
      scale(log(FlockSize)) +
      scale(log(BodyMass)) +
      scale(sin(rad)) + scale(cos(rad)) +
      # scale(Day)+
      scale(Temp) +
      scale(parks_percent_change_from_baseline) +
      (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality),
    data = ss[Nsp > 9], REML = FALSE
    )
    est_mg03c = est_out(mg03c, "03c) (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality); >9/specie")

    # export
      save(file = here::here("Data/Fig_S2_estimates.Rdata"), 
      est_m1a, est_m1b, est_m1c, est_m1d, est_m1e, est_m2a, est_m2b, est_m2c, est_m3a, est_m3b, est_m3c, 
      est_m01a, est_m01b, est_m01c, est_m02a, est_m02b, est_m02c, est_m03a, est_m03b, est_m03c, 
      est_mg01a, est_mg01b, est_mg01c, est_mg02a, est_mg02b, est_mg02c, est_mg03a, est_mg03b, est_mg03c)
     
     # prepare plot for Period
     xs = rbind(est_m1a, est_m1b, est_m1c, est_m1d, est_m1e, est_m2a, est_m2b, est_m2c, est_m3a, est_m3b, est_m3c)
     #gsub("scale\\(Covid\\)", "Period", "(scale(Covid)|Country)")
     xs[, model := gsub("scale\\(Covid\\)", "Period", xs$model)]
     xs[, model := gsub("Year", "year", xs$model)]
     xs[, model := gsub("Species", "species", xs$model)]
     xs[, model := gsub("sp_day_year", "species within day & year", xs$model)]
     xs[, model := gsub("IDLocality", "site", xs$model)]
     xs[, model := gsub("sp_loc", "species within site", xs$model)]

     gs2_ =
       ggplot(xs[predictor == "scale(Covid)"], aes(y = model, x = estimate, col = model)) +
       geom_vline(xintercept = 0, col = "grey30", lty = 3) +
       geom_errorbar(aes(xmin = lwr, xmax = upr, col = model), width = 0, position = position_dodge(width = 0.01)) +
       # ggtitle ("Sim based")+
       geom_point(position = position_dodge(width = 0.01)) +
       # scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
       # scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
       scale_color_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
       scale_fill_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
       scale_y_discrete(limits = rev) +
       coord_fixed(ratio = 0.05, xlim = c(-0.23, 0.15)) +
       # scale_shape(guide = guide_legend(reverse = TRUE)) +
       # scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
       labs(y = NULL, x = "Period (before vs during shutdowns)", tag = "a)") + # title = "a)",Effect of Period (before/during shutdown)
       # ylim(c(0,100))+
       # coord_flip()+
       theme_bw() +
       theme(
         legend.position = "none",
         #plot.title.position = "plot",
         plot.subtitle = element_text(size = 7),
         plot.title = element_text(size = 7),
         plot.tag = element_text(size = 7),
         legend.title = element_text(size = 7),
         legend.text = element_text(size = 6),
         ## legend.spacing.y = unit(0.1, 'cm'),
         legend.key.height = unit(0.5, "line"),
         # plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit =  "pt"),
         panel.grid = element_blank(),
         panel.border = element_blank(),
         panel.background = element_blank(),
         axis.line = element_line(colour = ax_lines, size = 0.25),
         axis.line.y = element_blank(),
         axis.ticks.y = element_blank(),
         axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
         axis.ticks.length = unit(1, "pt"),
         axis.text.x = element_text(colour = "black", size = 6),
         axis.text.y = element_text(colour = "black", size = 7),
         axis.title = element_text(size = 7)
       )
     #gs2
     # ggsave(here::here('Outputs/Figure_Sy.png'),g, width = 30, height =5, units = 'cm')
     # prepare plot for Stringency
     xs0 = rbind(est_m01a, est_m01b, est_m01c, est_m02a, est_m02b, est_m02c, est_m03a, est_m03b, est_m03c)
     xs0[, model := gsub("scale\\(StringencyIndex\\)", "Stringency Index", model)]
     xs0[, model := gsub("Year", "year", model)]
     xs0[, model := gsub("Species", "species", model)]
     xs0[, model := gsub("sp_day_year", "species within day & year", model)]
     xs0[, model := gsub("IDLocality", "site", model)]
     xs0[, model := gsub("sp_loc", "species within site", model)]
     g0_ =
       ggplot(xs0[predictor == "scale(StringencyIndex)"], aes(y = model, x = estimate, col = model)) +
       geom_vline(xintercept = 0, col = "grey30", lty = 3) +
       geom_errorbar(aes(xmin = lwr, xmax = upr, col = model), width = 0, position = position_dodge(width = 0.01)) +
       # ggtitle ("Sim based")+
       geom_point(position = position_dodge(width = 0.01)) +
       # scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
       # scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
       scale_color_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
       scale_fill_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
       scale_y_discrete(limits = rev) +
       coord_fixed(ratio = 0.05, xlim = c(-0.23, 0.15)) +
       # scale_shape(guide = guide_legend(reverse = TRUE)) +
       # scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
       labs(y = NULL, x = "Stringency index", tag = "b)") + # title = "b) Effect of ") +
       # ylim(c(0,100))+
       # coord_flip()+
       theme_bw() +
       theme(
         legend.position = "none",
         plot.subtitle = element_text(size = 7),
         plot.title = element_text(size = 7),
         plot.tag = element_text(size = 7),
         legend.title = element_text(size = 7),
         legend.text = element_text(size = 6),
         ## legend.spacing.y = unit(0.1, 'cm'),
         legend.key.height = unit(0.5, "line"),
         # plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit =  "pt"),
         panel.grid = element_blank(),
         panel.border = element_blank(),
         panel.background = element_blank(),
         axis.line = element_line(colour = ax_lines, size = 0.25),
         axis.line.y = element_blank(),
         axis.ticks.y = element_blank(),
         axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
         axis.ticks.length = unit(1, "pt"),
         axis.text.x = element_text(colour = "black", size = 6),
         axis.text.y = element_text(colour = "black", size = 7),
         axis.title = element_text(size = 7)
       )
     #g0
     # ggsave(here::here('Outputs/Figure_Sz.png'),g0, width = 30, height =5, units = 'cm')
     # prepare plot for Google
     xg0 = rbind(est_mg01a, est_mg01b, est_mg01c, est_mg02a, est_mg02b, est_mg02c, est_mg03a, est_mg03b, est_mg03c)
     xg0[, model := gsub("parks_percent_change_from_baseline", "Google Mobility", model)]
     xg0[, model := gsub("Year", "year", model)]
     xg0[, model := gsub("Species", "species", model)]
     xg0[, model := gsub("sp_day_year", "species within day & year", model)]
     xg0[, model := gsub("IDLocality", "site", model)]
     xg0[, model := gsub("sp_loc", "species within site", model)]
     gg0_ =
       ggplot(xg0[predictor == "scale(parks_percent_change_from_baseline)"], aes(y = model, x = estimate, col = model)) +
       geom_vline(xintercept = 0, col = "grey30", lty = 3) +
       geom_errorbar(aes(xmin = lwr, xmax = upr, col = model), width = 0, position = position_dodge(width = 0.01)) +
       # ggtitle ("Sim based")+
       geom_point(position = position_dodge(width = 0.01)) +
       # scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
       # scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
       scale_color_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
       scale_fill_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
       scale_y_discrete(limits = rev) +
       coord_fixed(ratio = 0.05, xlim = c(-0.23, 0.15)) +
       # scale_shape(guide = guide_legend(reverse = TRUE)) +
       # scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
       labs(y = NULL, x = "Google Mobility\n[Standardized effect sizes on\nflight initiation distances]", tag = "c)") + # title = "b) Effect of ") +
       # ylim(c(0,100))+
       # coord_flip()+
       theme_bw() +
       theme(
         legend.position = "none",
         plot.subtitle = element_text(size = 7),
         plot.title = element_text(size = 7),
         plot.tag = element_text(size = 7),
         legend.title = element_text(size = 7),
         legend.text = element_text(size = 6),
         ## legend.spacing.y = unit(0.1, 'cm'),
         legend.key.height = unit(0.5, "line"),
         # plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit =  "pt"),
         panel.grid = element_blank(),
         panel.border = element_blank(),
         panel.background = element_blank(),
         axis.line = element_line(colour = ax_lines, size = 0.25),
         axis.line.y = element_blank(),
         axis.ticks.y = element_blank(),
         axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
         axis.ticks.length = unit(1, "pt"),
         axis.text.x = element_text(colour = "black", size = 6),
         axis.text.y = element_text(colour = "black", size = 7),
         axis.title = element_text(size = 7)
       )
     #gg0
     # ggsave(here::here('Outputs/Figure_Sz.png'),g0, width = 30, height =5, units = 'cm')
     # combine
     grid.draw(rbind(ggplotGrob(gs2_), ggplotGrob(g0_), ggplotGrob(gg0_)))

     if(save_plot==TRUE){
     ggsave(here::here("Outputs/Fig_S2_rev_v6.png"), rbind(ggplotGrob(gs2_), ggplotGrob(g0_), ggplotGrob(gg0_)), width = 30, height = 15, units = "cm")
     }

Supplementary TABLES

Table S1 | Alternative models on escape distance given Period

    m1a_ = m_out(name = "Table S1 - 1a", dep = "Escape distance", model = m1a, nsim = 5000)
    m1b_ = m_out(name = "Table S1 - 1b", dep = "Escape distance", model = m1b, nsim = 5000)
    m1c_ = m_out(name = "Table S1 - 1c", dep = "Escape distance", model = m1c, nsim = 5000)
    m1d_ = m_out(name = "Table S1 - 1d", dep = "Escape distance", model = m1d, nsim = 5000)
    m1e_ = m_out(name = "Table S1 - 1e", dep = "Escape distance", model = m1e, nsim = 5000)

    m2a_ = m_out(name = "Table S1 - 2a", dep = "Escape distance", model = m2a, nsim = 5000)
    m2b_ = m_out(name = "Table S1 - 2b", dep = "Escape distance", model = m2b, nsim = 5000)
    m2c_ = m_out(name = "Table S1 - 2c", dep = "Escape distance", model = m2c, nsim = 5000)
    m3a_ = m_out(name = "Table S1 - 3a", dep = "Escape distance", model = m3a, nsim = 5000)
    m3b_ = m_out(name = "Table S1 - 3b", dep = "Escape distance", model = m3b, nsim = 5000)
    m3c_ = m_out(name = "Table S1 - 3c", dep = "Escape distance", model = m3c, nsim = 5000)
     
    out1 = rbind(m1a_, m1b_, m1c_, m1d_, m1e_, m2a_, m2b_, m2c_, m3a_, m3b_, m3c_, fill = TRUE)
    out1[is.na(out1)] = ""
    out1[, effect := gsub("Covid", "Period", effect)]
    out1[, effect := gsub("Year", "year", effect)]
    out1[, effect := gsub("Species", "species", effect)]
    out1[, effect := gsub("sp_day_year", "species within day & year", effect)]
    out1[, effect := gsub("IDLocality", "site", effect)]
    out1[, effect := gsub("sp_loc", "species within site", effect)]
    out1$R2_mar=out1$R2_con=NULL
    
    fwrite(file = here::here("Outputs/Table_S1_rev.csv"), out1)
    
    out1 %>%
       kbl() %>%
       kable_paper("hover", full_width = F)
model response error_structure N type effect estimate_r lwr_r upr_r
Table S1 - 1a Escape distance Gaussian 6369 fixed (Intercept) 0.082 -0.128 0.288
fixed scale(log(SD)) 0.477 0.454 0.501
fixed scale(log(FlockSize)) -0.008 -0.027 0.011
fixed scale(log(BodyMass)) 0.046 -0.02 0.114
fixed scale(sin(rad)) 0.022 -0.005 0.049
fixed scale(cos(rad)) 0.029 0.006 0.051
fixed scale(Temp) 0 -0.029 0.028
fixed scale(Period) -0.033 -0.191 0.13
random % species within day & year (Intercept) 8% 6% 9%
random % species within site (Intercept) 5% 4% 5%
random % site (Intercept) 1% 0% 5%
random % site scale(Period) 1% 0% 5%
random % species (Intercept) 9% 8% 8%
random % genus (Intercept) 6% 6% 6%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 1% -3% 5%
random % Country scale(Period) 1% -3% 5%
random % year (Intercept) 3% 1% 4%
random % Residual 65% 51% 76%
Table S1 - 1b Escape distance Gaussian 6369 fixed (Intercept) 0.069 -0.154 0.291
fixed scale(log(SD)) 0.478 0.455 0.503
fixed scale(log(FlockSize)) -0.007 -0.026 0.012
fixed scale(log(BodyMass)) 0.046 -0.019 0.112
fixed scale(sin(rad)) 0.02 -0.009 0.049
fixed scale(cos(rad)) 0.027 0.005 0.051
fixed scale(Temp) -0.001 -0.029 0.028
fixed scale(Period) -0.047 -0.209 0.111
random % species within day & year (Intercept) 12% 7% 14%
random % species within site (Intercept) 0% -1% 0%
random % species within site scale(Period) 0% -1% 0%
random % site (Intercept) 2% 0% 5%
random % site scale(Period) 2% 0% 5%
random % species (Intercept) 0% -1% 8%
random % species scale(Period) 0% -1% 8%
random % genus (Intercept) 0% 0% 5%
random % genus scale(Period) 0% 0% 5%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 1% -4% 6%
random % Country scale(Period) 1% -4% 6%
random % year (Intercept) 3% 2% 4%
random % Residual 79% 42% 94%
Table S1 - 1c Escape distance Gaussian 6369 fixed (Intercept) 0.013 -0.161 0.175
fixed scale(log(SD)) 0.494 0.469 0.518
fixed scale(log(FlockSize)) -0.018 -0.037 0
fixed scale(log(BodyMass)) 0.025 -0.006 0.054
fixed scale(sin(rad)) 0.014 -0.014 0.043
fixed scale(cos(rad)) 0.027 0.003 0.051
fixed scale(Temp) -0.004 -0.034 0.027
fixed scale(Period) -0.041 -0.212 0.132
random % species within day & year (Intercept) 12% 7% 13%
random % species within site (Intercept) 1% 0% 12%
random % species within site scale(Period) 1% 0% 12%
random % site (Intercept) 1% 1% 6%
random % site scale(Period) 1% 1% 6%
random % species scale(Period) 1% 1% 1%
random % genus scale(Period) 0% 0% 0%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 1% -2% 4%
random % Country scale(Period) 1% -2% 4%
random % year (Intercept) 3% 2% 4%
random % Residual 77% 46% 87%
Table S1 - 1d Escape distance Gaussian 6369 fixed (Intercept) 0.083 -0.124 0.286
fixed scale(log(SD)) 0.476 0.452 0.501
fixed scale(log(FlockSize)) -0.008 -0.026 0.011
fixed scale(log(BodyMass)) 0.046 -0.018 0.114
fixed scale(sin(rad)) 0.022 -0.007 0.049
fixed scale(cos(rad)) 0.028 0.006 0.051
fixed scale(Temp) 0 -0.029 0.029
fixed scale(Period) -0.031 -0.199 0.129
random % species within day & year (Intercept) 8% 6% 9%
random % species within site (Intercept) 0% 0% 4%
random % species within site scale(Period) 0% 0% 4%
random % site (Intercept) 1% 0% 5%
random % site scale(Period) 1% 0% 5%
random % species (Intercept) 10% 9% 10%
random % genus (Intercept) 7% 6% 6%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 1% -3% 5%
random % Country scale(Period) 1% -3% 5%
random % year (Intercept) 3% 2% 4%
random % Residual 66% 48% 78%
Table S1 - 1e Escape distance Gaussian 6369 fixed (Intercept) 0.077 -0.127 0.283
fixed scale(log(SD)) 0.478 0.454 0.501
fixed scale(log(FlockSize)) -0.008 -0.027 0.011
fixed scale(log(BodyMass)) 0.048 -0.021 0.114
fixed scale(sin(rad)) 0.021 -0.006 0.048
fixed scale(cos(rad)) 0.024 0.002 0.047
fixed scale(Temp) -0.006 -0.033 0.022
fixed scale(Period) -0.029 -0.204 0.142
random % species within day & year (Intercept) 7% 6% 8%
random % species within site (Intercept) 5% 4% 5%
random % site (Intercept) 6% 6% 6%
random % species (Intercept) 8% 8% 9%
random % genus (Intercept) 6% 5% 7%
random % weekday (Intercept) 0% 0% 1%
random % Country (Intercept) 1% -3% 5%
random % Country scale(Period) 1% -3% 5%
random % year (Intercept) 2% 1% 4%
random % Residual 62% 52% 71%
Table S1 - 2a Escape distance Gaussian 5261 fixed (Intercept) 0.161 -0.185 0.509
fixed scale(log(SD)) 0.463 0.437 0.49
fixed scale(log(FlockSize)) -0.01 -0.031 0.011
fixed scale(log(BodyMass)) -0.016 -0.123 0.095
fixed scale(sin(rad)) 0.027 -0.002 0.057
fixed scale(cos(rad)) 0.036 0.011 0.061
fixed scale(Temp) -0.005 -0.039 0.026
fixed scale(Period) -0.035 -0.204 0.133
random % species within day & year (Intercept) 7% 4% 9%
random % species within site (Intercept) 3% 2% 4%
random % site (Intercept) 1% -1% 2%
random % site scale(Period) 1% -1% 2%
random % species (Intercept) 10% 8% 9%
random % genus (Intercept) 5% 4% 4%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 3% -5% 20%
random % Country scale(Period) 3% -5% 20%
random % year (Intercept) 3% 2% 3%
random % Residual 64% 34% 84%
Table S1 - 2b Escape distance Gaussian 5261 fixed (Intercept) 0.164 -0.171 0.507
fixed scale(log(SD)) 0.462 0.435 0.489
fixed scale(log(FlockSize)) -0.01 -0.029 0.01
fixed scale(log(BodyMass)) -0.011 -0.122 0.096
fixed scale(sin(rad)) 0.028 -0.003 0.057
fixed scale(cos(rad)) 0.036 0.011 0.061
fixed scale(Temp) -0.005 -0.037 0.028
fixed scale(Period) -0.036 -0.204 0.146
random % species within day & year (Intercept) 9% 4% 11%
random % species within site (Intercept) 0% 0% 2%
random % species within site scale(Period) 0% 0% 2%
random % site (Intercept) 1% -1% 2%
random % site scale(Period) 1% -1% 2%
random % species (Intercept) 0% -1% 5%
random % species scale(Period) 0% -1% 5%
random % genus (Intercept) 0% 0% 6%
random % genus scale(Period) 0% 0% 6%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 3% -6% 17%
random % Country scale(Period) 3% -6% 17%
random % year (Intercept) 4% 2% 3%
random % Residual 78% 30% 100%
Table S1 - 2c Escape distance Gaussian 5261 fixed (Intercept) 0.153 -0.197 0.513
fixed scale(log(SD)) 0.464 0.437 0.492
fixed scale(log(FlockSize)) -0.01 -0.03 0.011
fixed scale(log(BodyMass)) -0.01 -0.114 0.101
fixed scale(sin(rad)) 0.026 -0.002 0.055
fixed scale(cos(rad)) 0.032 0.007 0.057
fixed scale(Temp) -0.012 -0.044 0.02
fixed scale(Period) -0.031 -0.213 0.147
random % species within day & year (Intercept) 7% 5% 9%
random % species within site (Intercept) 3% 2% 4%
random % site (Intercept) 4% 3% 4%
random % species (Intercept) 8% 7% 7%
random % genus (Intercept) 6% 5% 6%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 2% -6% 17%
random % Country scale(Period) 2% -6% 17%
random % year (Intercept) 3% 2% 4%
random % Residual 63% 39% 81%
Table S1 - 3a Escape distance Gaussian 5107 fixed (Intercept) 0.137 -0.205 0.493
fixed scale(log(SD)) 0.462 0.435 0.488
fixed scale(log(FlockSize)) -0.01 -0.031 0.01
fixed scale(log(BodyMass)) 0.003 -0.104 0.114
fixed scale(sin(rad)) 0.028 -0.003 0.058
fixed scale(cos(rad)) 0.034 0.01 0.06
fixed scale(Temp) -0.006 -0.04 0.028
fixed scale(Period) -0.032 -0.207 0.142
random % species within day & year (Intercept) 7% 4% 9%
random % species within site (Intercept) 3% 2% 4%
random % site (Intercept) 1% -1% 2%
random % site scale(Period) 1% -1% 2%
random % species (Intercept) 21% 16% 18%
random % genus (Intercept) 0% 0% 0%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 2% -5% 19%
random % Country scale(Period) 2% -5% 19%
random % year (Intercept) 3% 2% 3%
random % Residual 60% 32% 78%
Table S1 - 3b Escape distance Gaussian 5107 fixed (Intercept) 0.139 -0.21 0.484
fixed scale(log(SD)) 0.461 0.433 0.487
fixed scale(log(FlockSize)) -0.01 -0.031 0.01
fixed scale(log(BodyMass)) 0.011 -0.102 0.126
fixed scale(sin(rad)) 0.027 -0.002 0.057
fixed scale(cos(rad)) 0.034 0.009 0.06
fixed scale(Temp) -0.006 -0.04 0.028
fixed scale(Period) -0.03 -0.201 0.142
random % species within day & year (Intercept) 9% 3% 11%
random % species within site (Intercept) 0% 0% 2%
random % species within site scale(Period) 0% 0% 2%
random % site (Intercept) 1% -1% 2%
random % site scale(Period) 1% -1% 2%
random % species (Intercept) 0% -1% 9%
random % species scale(Period) 0% -1% 9%
random % genus (Intercept) 0% 0% 3%
random % genus scale(Period) 0% 0% 3%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 3% -5% 16%
random % Country scale(Period) 3% -5% 16%
random % year (Intercept) 4% 2% 3%
random % Residual 78% 28% 100%
Table S1 - 3c Escape distance Gaussian 5107 fixed (Intercept) 0.127 -0.238 0.486
fixed scale(log(SD)) 0.463 0.437 0.49
fixed scale(log(FlockSize)) -0.01 -0.031 0.01
fixed scale(log(BodyMass)) 0.005 -0.101 0.113
fixed scale(sin(rad)) 0.026 -0.003 0.056
fixed scale(cos(rad)) 0.029 0.004 0.053
fixed scale(Temp) -0.014 -0.046 0.018
fixed scale(Period) -0.027 -0.205 0.15
random % species within day & year (Intercept) 7% 5% 9%
random % species within site (Intercept) 3% 2% 4%
random % site (Intercept) 5% 3% 5%
random % species (Intercept) 15% 11% 15%
random % genus (Intercept) 0% 0% 0%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 2% -6% 18%
random % Country scale(Period) 2% -6% 18%
random % year (Intercept) 3% 2% 4%
random % Residual 62% 38% 78%

Continuous variables were z-transformed (scale). Note that (1a) model is the one reported in the main text.

Table S2 | Alternative models on escape distance given Stringency

    m01a_ = m_out(name = "Table S2 - 1a", dep = "Escape distance", model = m01a, nsim = 5000)
    m01b_ = m_out(name = "Table S2 - 1b", dep = "Escape distance", model = m01b, nsim = 5000)
    m01c_ = m_out(name = "Table S2 - 1c", dep = "Escape distance", model = m01c, nsim = 5000)

    m02a_ = m_out(name = "Table S2 - 2a", dep = "Escape distance", model = m02a, nsim = 5000)
    m02b_ = m_out(name = "Table S2 - 2b", dep = "Escape distance", model = m02b, nsim = 5000)
    m02c_ = m_out(name = "Table S2 - 2c", dep = "Escape distance", model = m02c, nsim = 5000)

    m03a_ = m_out(name = "Table S2 - 3a", dep = "Escape distance", model = m03a, nsim = 5000)
    m03b_ = m_out(name = "Table S2 - 3b", dep = "Escape distance", model = m03b, nsim = 5000)
    m03c_ = m_out(name = "Table S2 - 3c", dep = "Escape distance", model = m03c, nsim = 5000)

    out2 = rbind(m01a_, m01b_, m01c_, m02a_, m02b_, m02c_, m03a_, m03b_, m03c_, fill = TRUE)
    out2[is.na(out2)] = ""
    out2[, effect := gsub("StringencyIndex", "Stringency Index", effect)]
    out2[, effect := gsub("Year", "year", effect)]
    out2[, effect := gsub("Species", "species", effect)]
    out2[, effect := gsub("sp_day_year", "species within day & year", effect)]
    out2[, effect := gsub("IDLocality", "site", effect)]
    out2[, effect := gsub("sp_loc", "species within site", effect)]
    out2$R2_mar = out2$R2_con = NULL

    fwrite(file = here::here("Outputs/Table_S2_rev.csv"), out2)

    out2 %>%
      kbl() %>%
      kable_paper("hover", full_width = F)
model response error_structure N type effect estimate_r lwr_r upr_r
Table S2 - 1a Escape distance Gaussian 3676 fixed (Intercept) -0.075 -0.329 0.177
fixed scale(year) 0.019 -0.03 0.065
fixed scale(log(SD)) 0.492 0.46 0.523
fixed scale(log(FlockSize)) -0.004 -0.029 0.024
fixed scale(log(BodyMass)) 0.022 -0.046 0.091
fixed scale(sin(rad)) -0.02 -0.061 0.022
fixed scale(cos(rad)) -0.005 -0.04 0.029
fixed scale(Temp) -0.007 -0.048 0.033
fixed scale(Stringency Index) 0.014 -0.157 0.178
random % species within day & year (Intercept) 5% 3% 7%
random % species within site (Intercept) 5% 3% 6%
random % species (Intercept) 7% 6% 9%
random % site (Intercept) 8% 7% 9%
random % genus (Intercept) 4% 4% 4%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 2% -15% 16%
random % Country scale(Stringency Index) 2% -15% 16%
random % Residual 67% 45% 95%
Table S2 - 1b Escape distance Gaussian 3676 fixed (Intercept) -0.074 -0.329 0.181
fixed scale(year) 0.019 -0.029 0.066
fixed scale(log(SD)) 0.493 0.461 0.525
fixed scale(log(FlockSize)) -0.003 -0.029 0.023
fixed scale(log(BodyMass)) 0.018 -0.051 0.087
fixed scale(sin(rad)) -0.02 -0.06 0.021
fixed scale(cos(rad)) -0.004 -0.041 0.03
fixed scale(Temp) -0.008 -0.051 0.035
fixed scale(Stringency Index) 0.008 -0.161 0.17
random % species within day & year (Intercept) 5% 3% 7%
random % species within site (Intercept) 5% 3% 7%
random % species (Intercept) 7% 5% 8%
random % site (Intercept) 8% 6% 10%
random % genus (Intercept) 0% -1% 4%
random % genus scale(Stringency Index) 0% -1% 4%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 2% -16% 15%
random % Country scale(Stringency Index) 2% -16% 15%
random % Residual 70% 44% 100%
Table S2 - 1c Escape distance Gaussian 3676 fixed (Intercept) -0.109 -0.277 0.06
fixed scale(year) -0.003 -0.051 0.045
fixed scale(log(SD)) 0.55 0.518 0.583
fixed scale(log(FlockSize)) -0.023 -0.049 0.003
fixed scale(log(BodyMass)) -0.024 -0.053 0.005
fixed scale(sin(rad)) -0.039 -0.082 0.002
fixed scale(cos(rad)) 0.005 -0.032 0.041
fixed scale(Temp) 0.003 -0.039 0.045
fixed scale(Stringency Index) 0.019 -0.134 0.176
random % site (Intercept) 10% 10% 10%
random % Country (Intercept) 1% -7% 8%
random % Country scale(Stringency Index) 1% -7% 8%
random % Residual 88% 74% 103%
Table S2 - 2a Escape distance Gaussian 3573 fixed (Intercept) -0.153 -0.445 0.134
fixed scale(year) 0.022 -0.026 0.071
fixed scale(log(SD)) 0.483 0.451 0.515
fixed scale(log(FlockSize)) 0.002 -0.024 0.03
fixed scale(log(BodyMass)) -0.024 -0.099 0.053
fixed scale(sin(rad)) -0.02 -0.061 0.021
fixed scale(cos(rad)) -0.008 -0.044 0.028
fixed scale(Temp) -0.008 -0.049 0.034
fixed scale(Stringency Index) 0.021 -0.15 0.192
random % species within day & year (Intercept) 5% 3% 7%
random % species within site (Intercept) 4% 3% 7%
random % site (Intercept) 7% 5% 9%
random % species (Intercept) 11% 9% 13%
random % genus (Intercept) 2% 2% 2%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 3% -20% 20%
random % Country scale(Stringency Index) 3% -20% 20%
random % Residual 66% 39% 102%
Table S2 - 2b Escape distance Gaussian 3573 fixed (Intercept) -0.153 -0.439 0.135
fixed scale(year) 0.021 -0.028 0.07
fixed scale(log(SD)) 0.484 0.451 0.517
fixed scale(log(FlockSize)) 0.003 -0.024 0.031
fixed scale(log(BodyMass)) -0.03 -0.107 0.046
fixed scale(sin(rad)) -0.019 -0.061 0.022
fixed scale(cos(rad)) -0.006 -0.042 0.029
fixed scale(Temp) -0.008 -0.051 0.033
fixed scale(Stringency Index) 0.019 -0.152 0.185
random % species within day & year (Intercept) 5% 3% 8%
random % species within site (Intercept) 4% 3% 7%
random % site (Intercept) 7% 5% 9%
random % species (Intercept) 11% 8% 13%
random % genus (Intercept) 0% -1% 2%
random % genus scale(Stringency Index) 0% -1% 2%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 3% -21% 19%
random % Country scale(Stringency Index) 3% -21% 19%
random % Residual 67% 38% 106%
Table S2 - 2c Escape distance Gaussian 3573 fixed (Intercept) -0.107 -0.288 0.07
fixed scale(year) -0.001 -0.049 0.048
fixed scale(log(SD)) 0.543 0.51 0.575
fixed scale(log(FlockSize)) -0.019 -0.045 0.009
fixed scale(log(BodyMass)) -0.027 -0.056 0.003
fixed scale(sin(rad)) -0.039 -0.08 0.003
fixed scale(cos(rad)) 0.002 -0.036 0.04
fixed scale(Temp) 0.002 -0.04 0.046
fixed scale(Stringency Index) 0.024 -0.14 0.193
random % site (Intercept) 10% 10% 10%
random % Country (Intercept) 1% -7% 8%
random % Country scale(Stringency Index) 1% -7% 8%
random % Residual 88% 73% 105%
Table S2 - 3a Escape distance Gaussian 3425 fixed (Intercept) -0.134 -0.442 0.179
fixed scale(year) 0.027 -0.023 0.076
fixed scale(log(SD)) 0.481 0.448 0.515
fixed scale(log(FlockSize)) 0.002 -0.024 0.029
fixed scale(log(BodyMass)) 0.024 -0.067 0.11
fixed scale(sin(rad)) -0.011 -0.054 0.032
fixed scale(cos(rad)) -0.009 -0.045 0.028
fixed scale(Temp) -0.011 -0.054 0.032
fixed scale(Stringency Index) 0.018 -0.152 0.188
random % species within day & year (Intercept) 5% 3% 7%
random % species within site (Intercept) 5% 3% 7%
random % site (Intercept) 8% 5% 10%
random % species (Intercept) 9% 7% 10%
random % genus (Intercept) 3% 3% 3%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 2% -20% 20%
random % Country scale(Stringency Index) 2% -20% 20%
random % Residual 67% 39% 104%
Table S2 - 3b Escape distance Gaussian 3425 fixed (Intercept) -0.138 -0.449 0.171
fixed scale(year) 0.026 -0.025 0.076
fixed scale(log(SD)) 0.482 0.45 0.514
fixed scale(log(FlockSize)) 0.001 -0.025 0.029
fixed scale(log(BodyMass)) 0.02 -0.07 0.106
fixed scale(sin(rad)) -0.012 -0.054 0.031
fixed scale(cos(rad)) -0.009 -0.046 0.029
fixed scale(Temp) -0.01 -0.053 0.032
fixed scale(Stringency Index) 0.015 -0.157 0.194
random % species within day & year (Intercept) 5% 3% 7%
random % species within site (Intercept) 4% 3% 7%
random % site (Intercept) 8% 5% 11%
random % species (Intercept) 10% 7% 12%
random % genus (Intercept) 0% 0% 2%
random % genus scale(Stringency Index) 0% 0% 2%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 3% -23% 21%
random % Country scale(Stringency Index) 3% -23% 21%
random % Residual 67% 36% 111%
Table S2 - 3c Escape distance Gaussian 3425 fixed (Intercept) -0.106 -0.292 0.071
fixed scale(year) 0.005 -0.044 0.055
fixed scale(log(SD)) 0.541 0.508 0.575
fixed scale(log(FlockSize)) -0.02 -0.046 0.007
fixed scale(log(BodyMass)) -0.019 -0.049 0.011
fixed scale(sin(rad)) -0.031 -0.074 0.012
fixed scale(cos(rad)) -0.001 -0.039 0.039
fixed scale(Temp) -0.002 -0.045 0.042
fixed scale(Stringency Index) 0.022 -0.145 0.18
random % site (Intercept) 11% 11% 11%
random % Country (Intercept) 1% -8% 8%
random % Country scale(Stringency Index) 1% -8% 8%
random % Residual 88% 73% 105%

Continuous variables were z-transformed (scale). Note that (1a) model is the one reported in the main text.

Table S3 | Alternative models on escape distance given Google Mobility

    # google
     mg01a_ = m_out(name = "Table S3 - 1a", dep = "Escape distance", model = mg01a, nsim = 5000)
     mg01b_ = m_out(name = "Table S3 - 1b", dep = "Escape distance", model = mg01b, nsim = 5000)
     mg01c_ = m_out(name = "Table S3 - 1c", dep = "Escape distance", model = mg01c, nsim = 5000)
## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.
     mg02a_ = m_out(name = "Table S3 - 2a", dep = "Escape distance", model = mg02a, nsim = 5000)
     mg02b_ = m_out(name = "Table S3 - 2b", dep = "Escape distance", model = mg02b, nsim = 5000)
     mg02c_ = m_out(name = "Table S3 - 2c", dep = "Escape distance", model = mg02c, nsim = 5000)
## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.
     mg03a_ = m_out(name = "Table S3 - 3a", dep = "Escape distance", model = mg03a, nsim = 5000)
     mg03b_ = m_out(name = "Table S3 - 3b", dep = "Escape distance", model = mg03b, nsim = 5000)
     mg03c_ = m_out(name = "Table S3 - 3c", dep = "Escape distancey", model = mg03c, nsim = 5000)
## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.
     out3 = rbind(mg01a_, mg01b_, mg01c_, mg02a_, mg02b_, mg02c_, mg03a_, mg03b_, mg03c_, fill = TRUE)
     out3[is.na(out3)] = ""
       out3[, effect := gsub("parks_percent_change_from_baseline", "Google Mobility", effect)]
       out3[, effect := gsub("Year", "year", effect)]
       out3[, effect := gsub("Species", "species", effect)]
       out3[, effect := gsub("sp_day_year", "species within day & year", effect)]
       out3[, effect := gsub("IDLocality", "site", effect)]
       out3[, effect := gsub("sp_loc", "species within site", effect)]
       out3$R2_mar = out3$R2_con = NULL
     fwrite(file = here::here("Outputs/Table_S3_rev.csv"), out3)
     
      out3 %>%
          kbl() %>%
          kable_paper("hover", full_width = F)
model response error_structure N type effect estimate_r lwr_r upr_r
Table S3 - 1a Escape distance Gaussian 3399 fixed (Intercept) -0.015 -0.169 0.141
fixed scale(year) 0.027 -0.017 0.071
fixed scale(log(SD)) 0.487 0.455 0.52
fixed scale(log(FlockSize)) 0.002 -0.024 0.029
fixed scale(log(BodyMass)) 0.026 -0.061 0.113
fixed scale(sin(rad)) 0.001 -0.04 0.044
fixed scale(cos(rad)) 0.008 -0.027 0.042
fixed scale(Temp) 0.022 -0.017 0.06
fixed scale(Google Mobility) -0.108 -0.232 0.016
random % species within day & year (Intercept) 6% 5% 6%
random % species within site (Intercept) 4% 4% 5%
random % site (Intercept) 9% 9% 9%
random % species (Intercept) 7% 6% 9%
random % genus (Intercept) 3% 3% 4%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 0% -6% 5%
random % Country scale(Google Mobility) 0% -6% 5%
random % Residual 70% 60% 82%
Table S3 - 1b Escape distance Gaussian 3399 fixed (Intercept) -0.013 -0.165 0.144
fixed scale(year) 0.026 -0.017 0.07
fixed scale(log(SD)) 0.487 0.456 0.522
fixed scale(log(FlockSize)) 0.002 -0.026 0.03
fixed scale(log(BodyMass)) 0.032 -0.051 0.119
fixed scale(sin(rad)) 0 -0.041 0.043
fixed scale(cos(rad)) 0.006 -0.028 0.042
fixed scale(Temp) 0.023 -0.017 0.064
fixed scale(Google Mobility) -0.111 -0.237 0.016
random % species within day & year (Intercept) 6% 5% 7%
random % species within site (Intercept) 5% 4% 5%
random % site (Intercept) 9% 9% 9%
random % species (Intercept) 9% 8% 9%
random % genus (Intercept) 0% -1% 4%
random % genus scale(Google Mobility) 0% -1% 4%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 0% -7% 5%
random % Country scale(Google Mobility) 0% -7% 5%
random % Residual 71% 56% 86%
Table S3 - 1c Escape distance Gaussian 3399 fixed (Intercept) -0.061 -0.15 0.024
fixed scale(year) 0.022 -0.017 0.062
fixed scale(log(SD)) 0.549 0.517 0.582
fixed scale(log(FlockSize)) -0.018 -0.046 0.008
fixed scale(log(BodyMass)) -0.019 -0.048 0.011
fixed scale(sin(rad)) -0.024 -0.065 0.019
fixed scale(cos(rad)) 0.015 -0.021 0.052
fixed scale(Temp) 0.03 -0.009 0.071
fixed scale(Google Mobility) -0.092 -0.206 0.026
random % site (Intercept) 11% 9% 12%
random % Country (Intercept) 1% -1% 3%
random % Country scale(Google Mobility) 1% -1% 3%
random % Residual 88% 82% 92%
Table S3 - 2a Escape distance Gaussian 3399 fixed (Intercept) -0.017 -0.17 0.138
fixed scale(year) 0.027 -0.017 0.07
fixed scale(log(SD)) 0.487 0.454 0.52
fixed scale(log(FlockSize)) 0.002 -0.024 0.03
fixed scale(log(BodyMass)) 0.026 -0.062 0.113
fixed scale(sin(rad)) 0.001 -0.04 0.043
fixed scale(cos(rad)) 0.007 -0.028 0.042
fixed scale(Temp) 0.022 -0.017 0.062
fixed scale(Google Mobility) -0.106 -0.234 0.018
random % species within day & year (Intercept) 6% 5% 6%
random % species within site (Intercept) 4% 4% 5%
random % site (Intercept) 9% 9% 9%
random % species (Intercept) 7% 6% 9%
random % genus (Intercept) 3% 3% 4%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 0% -6% 5%
random % Country scale(Google Mobility) 0% -6% 5%
random % Residual 70% 60% 83%
Table S3 - 2b Escape distance Gaussian 3399 fixed (Intercept) -0.014 -0.167 0.14
fixed scale(year) 0.027 -0.016 0.07
fixed scale(log(SD)) 0.488 0.455 0.521
fixed scale(log(FlockSize)) 0.002 -0.025 0.029
fixed scale(log(BodyMass)) 0.032 -0.053 0.115
fixed scale(sin(rad)) 0 -0.041 0.042
fixed scale(cos(rad)) 0.007 -0.029 0.041
fixed scale(Temp) 0.023 -0.018 0.063
fixed scale(Google Mobility) -0.108 -0.235 0.021
random % species within day & year (Intercept) 6% 5% 7%
random % species within site (Intercept) 5% 4% 5%
random % site (Intercept) 9% 9% 9%
random % species (Intercept) 9% 8% 9%
random % genus (Intercept) 0% -1% 4%
random % genus scale(Google Mobility) 0% -1% 4%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 0% -6% 5%
random % Country scale(Google Mobility) 0% -6% 5%
random % Residual 71% 56% 86%
Table S3 - 2c Escape distance Gaussian 3399 fixed (Intercept) -0.06 -0.145 0.03
fixed scale(year) 0.022 -0.018 0.063
fixed scale(log(SD)) 0.55 0.516 0.583
fixed scale(log(FlockSize)) -0.019 -0.046 0.008
fixed scale(log(BodyMass)) -0.018 -0.05 0.01
fixed scale(sin(rad)) -0.024 -0.066 0.018
fixed scale(cos(rad)) 0.015 -0.021 0.051
fixed scale(Temp) 0.03 -0.008 0.069
fixed scale(Google Mobility) -0.091 -0.206 0.023
random % site (Intercept) 11% 10% 12%
random % Country (Intercept) 1% -1% 3%
random % Country scale(Google Mobility) 1% -1% 3%
random % Residual 88% 82% 92%
Table S3 - 3a Escape distance Gaussian 3399 fixed (Intercept) -0.015 -0.169 0.141
fixed scale(year) 0.027 -0.016 0.071
fixed scale(log(SD)) 0.487 0.455 0.521
fixed scale(log(FlockSize)) 0.002 -0.024 0.029
fixed scale(log(BodyMass)) 0.026 -0.058 0.112
fixed scale(sin(rad)) 0.001 -0.042 0.042
fixed scale(cos(rad)) 0.008 -0.028 0.042
fixed scale(Temp) 0.022 -0.017 0.063
fixed scale(Google Mobility) -0.107 -0.235 0.02
random % species within day & year (Intercept) 6% 5% 6%
random % species within site (Intercept) 4% 4% 5%
random % site (Intercept) 9% 8% 9%
random % species (Intercept) 7% 6% 9%
random % genus (Intercept) 3% 3% 4%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 0% -6% 5%
random % Country scale(Google Mobility) 0% -6% 5%
random % Residual 70% 60% 82%
Table S3 - 3b Escape distance Gaussian 3399 fixed (Intercept) -0.011 -0.161 0.144
fixed scale(year) 0.027 -0.017 0.071
fixed scale(log(SD)) 0.488 0.454 0.521
fixed scale(log(FlockSize)) 0.002 -0.026 0.028
fixed scale(log(BodyMass)) 0.032 -0.052 0.116
fixed scale(sin(rad)) 0 -0.041 0.042
fixed scale(cos(rad)) 0.006 -0.028 0.042
fixed scale(Temp) 0.023 -0.016 0.063
fixed scale(Google Mobility) -0.109 -0.237 0.017
random % species within day & year (Intercept) 6% 5% 7%
random % species within site (Intercept) 5% 4% 5%
random % site (Intercept) 9% 9% 9%
random % species (Intercept) 9% 8% 9%
random % genus (Intercept) 0% -1% 4%
random % genus scale(Google Mobility) 0% -1% 4%
random % weekday (Intercept) 0% 0% 0%
random % Country (Intercept) 0% -6% 5%
random % Country scale(Google Mobility) 0% -6% 5%
random % Residual 71% 56% 86%
Table S3 - 3c Escape distancey Gaussian 3399 fixed (Intercept) -0.061 -0.149 0.026
fixed scale(year) 0.022 -0.017 0.064
fixed scale(log(SD)) 0.549 0.516 0.582
fixed scale(log(FlockSize)) -0.019 -0.046 0.008
fixed scale(log(BodyMass)) -0.019 -0.049 0.011
fixed scale(sin(rad)) -0.024 -0.064 0.018
fixed scale(cos(rad)) 0.015 -0.022 0.053
fixed scale(Temp) 0.03 -0.01 0.07
fixed scale(Google Mobility) -0.092 -0.205 0.024
random % site (Intercept) 11% 9% 12%
random % Country (Intercept) 1% -1% 3%
random % Country scale(Google Mobility) 1% -1% 3%
random % Residual 88% 82% 92%

Continuous variables were z-transformed (scale). Note that (1a) model is the one reported in the main text.

# TO DO modelAss

# END